Open quantum entanglement: a study of two atomic system in static patch of de Sitter space

Abstract

In this work, our prime objective is to study non-locality and long range effect of two body correlation using quantum entanglement from various information theoretic measure in the static patch of de Sitter space using a two body Open Quantum System (OQS). The OQS is described by a system of two entangled atoms, surrounded by a thermal bath, which is modelled by a massless probe scalar field. Firstly, we partially trace over the bath field and construct the Gorini Kossakowski Sudarshan Lindblad (GSKL) master equation, which describes the time evolution of the reduced subsystem density matrix. This GSKL master equation is characterized by two components, these are-Spin chain interaction Hamiltonian and the Lindbladian. To fix the form of both of them, we compute the Wightman functions for probe massless scalar field. Using this result alongwith the large time equilibrium behaviour we obtain the analytical solution for reduced density matrix. Further using this solution we evaluate various entanglement measures, namely Von-Neumann entropy, R$e'$nyi entropy, logarithmic negativity, entanglement of formation, concurrence and quantum discord for the two atomic subsystem on the static patch of De-Sitter space. Finally, we have studied violation of Bell-CHSH inequality, which is the key ingredient to study non-locality in primordial cosmology.

1 Introduction

The theory of closed quantum systems is a very popular topic and has already been firmly established. But in practical situations no quantum system can be ideally treated as closed and its interactions with the surroundings cannot be neglected, hence it becomes essential to develop a theoretical framework to treat these non adiabatic interactions and develop a proper understanding of the quantum mechanical system. The knowledge of complete time evolutionary dynamics of a quantum mechanical system requires incorporation of the details of the thermal environment, which paves the way towards the study of Open Quantum System [1, 43,44,45,46] (OQS), where the physical system weakly interacts with the environment. In general these interactions with the environment significantly controls the time evolution of the quantum mechanical system and induces the phenomenon of quantum dissipation and gives rise to many out of equilibrium related phenomenon [2, 3]. The dynamics of the reduced subsystem of OQS cannot be described using the unitary time evolution operators after integrating out the bath degrees of freedom from the theory. Correct time evolution of the system requires solving the effective master equation, which describes the non-unitary time evolution of the reduced density matrix of the system.

To deal with the subsystem in the context of OQS the entire combination of the system and its surroundings (thermal bath) is together treated as a closed quantum system and hence the evolution equations can be assumed to follow the unitary transformation rules. The prime assumption used in the present context is that the entire system environment combination forms a large closed quantum system. Therefore, its time evolution is governed by a unitary transformation generated by a global Hamiltonian which is made up of subsystem Hamiltonian, bath Hamiltonian and interaction Hamiltonian. Moreover the interaction between the system and the bath is assumed to depend on only the present moment, it carries no past memories at all. In short, the interaction is Markovian in nature. The interaction between the system and the surroundings is assumed to be weak which justifies the argument that the only effective change that can be seen over time occurs in the context of OQS. This assumption is generally useful in treating the evolution of the system when the OQS has sufficient time to relax to the equilibrium before being perturbed in presence of interaction between the system and thermal environment. However, in a special situation where the system has very fast or frequent perturbation in presence of system- thermal bath coupling, one needs to consider Non Markovian approximation. In this treatment additionally it has been assumed that the system is completely uncorrelated with the surroundings at initial time scale, provided the coupling between the subsystem and the environment is sufficiently weak in nature. The techniques developed in the context of OQS have proven very powerful in the context of quantum optics, statistical mechanics, information theory, thermodynamics, cosmology and biology.

On the other hand, quantum entanglement [4] is probably the most fascinating manifestation of quantum theory in which the beauty of quantum mechanics is truly realized. It is a physical phenomenon that occurs when pairs or groups of particles are generated, interact, or share spatial proximity in ways such that the quantum state of each particle cannot be described independently of the state of the others, even when the particles are separated by a large distance. Equivalently quantum entanglement of a state which is shared by two parties is necessary but not sufficient for that state to be non-local. Quantum entanglement plays a significant role in the context of quantum computation [5], information and teleportation theory [6, 7], quantum error correction [8, 9], etc. There are huge applications of quantum entanglement in various contexts. In interferometry the process of entanglement is used to achieve the Heisenberg limit [10]. In multi electron atomic system the electronic shells always consists of electrons in entangled state [11]. In the process of photosynthesis,entanglement is seen in the transfer of energy between light harvesting complexes and photosynthetic reaction centres [12]. In living organisms like bacteria entanglement has been observed between the organism and quantized light [13]. Study of various quantum information theoretic measures which quantifies the phenomenon of quantum entanglement from different types of OQS are important topics of research at present [14]. Few of them, namely Von Neumann entanglement entropy, R$e'$nyi entropy, quantum discord, concurrence and entanglement of formation are calculated in this paper to study the explicit role of quantum entanglement in the reduced subsystem between the two atoms.

Many authors have studied the physics of quantum fields in curved spacetime using one and two atomic system in the context of OQS [15,16,17]. A single detector system weakly interacting with a reservoir in quantized conformally coupled scalar field in De-Sitter space was investigated in the context of OQS in [18]. A similar study was done using two atoms (detectors) in [15] where the authors have studied the evolution of the subsystem only under the effect of the Lindbladian operator for different initial states. In that paper the authors have used only concurrence as the measure to quantify entanglement.In this present work we study the entanglement generation and dynamics, not only in the early and the late time scales but also at any arbitrary time scale. For that purpose we have treated the two atomic system as an open quantum system that interacts with the massless probe scalar field acting as the bath in the de Sitter curved background.After tracing out the bath degrees of freedom, we have analytically solved the GSKL master equation, which gives the non unitary evolution of any open quantum system in constant interaction with the bath. We have taken into account the contributions from both the prime components of the master equation, the effective Lamb Shift Hamiltonian and the Lindbladian operator, which provides a complete solution of the master equation and gives us a proper understanding about the complete time evolution of the reduced density matrix of the two atomic subsystem. Further using this result we have computed a number of information theoretic measures to quantify the quantum entanglement and study the entanglement dynamics and its dependence on various parameters of the system and background spacetime. Additionally, we have studied quantum non-locality by establishing Bell’s CHSH inequality violation in de Sitter space from the present OQS two atomic set up.

In Fig. 1, we have explicitly shown the Penrose diagram representing the static patch of the de Sitter space in which we have placed a two atomic open quantum system which is interacting with a thermal bath non-adiabatically. This diagram represents actually the causal patch of an observer sitting at the north pole, which is represented by $r=0$ is static coordinate in de Sitter space. Equivalently this can be described in global coordinate with $\theta =0$. Here the bifurcation Killing horizon for $\partial _{t}$ is represented by $r=\alpha$ where the parameter $\alpha =\sqrt{\frac{3}{\Lambda }}>0$ as the cosmological constant $\Lambda >0$ in de Sitter space. In this context the bifurcation sphere appears just at the middle of the Penrose diagram and in global coordinates it is described by $t=0$ time slice. However, the other three regions can also be filled by the static coordinate system, which is just like Schwarzschild black holes. In Fig. 2, we have presented a mnemonic chart of the computational scheme for the present two atomic OQS set up.

Fig. 1

Penrose diagram representing the static patch of the de Sitter space in which we have placed a two atomic open quantum system which is interacting with a thermal bath non-adiabatically. This diagram actually represents the causal patch of an observer sitting at the north pole, which is represented by $r=0$ is static coordinate in de Sitter space. Equivalently this can be described in global coordinate with $\theta =0$. Here the bifurcation Killing horizon for $\partial _{t}$ is represented by $r=\alpha$ where the parameter $\alpha =\sqrt{\frac{3}{\Lambda }}>0$ as the cosmological constant $\Lambda >0$ in de Sitter space. In this context the bifurcation sphere appears just at the middle of the Penrose diagram and in global coordinates it is described by $t=0$ time slice. However, the other three regions can also be filled by the static coordinate system, which is just like Schwarzschild black holes

Fig. 2

Mnemonic chart for the two atomic (two body) entangled OQS set up

The plan of this paper is as follows. In Sect. 2, we discuss the basics of the OQS. The objective of this section is to familiarize the reader with the prime components of open quantum systems. Further, in Sect. 4 and Sect. 5, we show the explicit construction of the two prime components of Gorini Kossakowski Sudarshan Lindblad (GSKL) master equation namely the effective Lamb Shift part of the Hamiltonian and the quantum dissipator operator or the Lindbladian. In Sect. 6, we explicitly calculate the analytical solution of the Gorini Kossakowski Sudarshan Lindblad (GSKL) master equation  [19, 20] for the case of two entangled atoms which mimics the role of Unruh-De-Witt detectors, which are minimally coupled to a probe massless scalar field placed in the thermal bath. In Sects. 7.1–7.5, we explicitly calculate various entanglement measures used in the context of quantum information theory now-a-days. To name some we calculate Von Neumann entanglement entropy [21], R$e'$nyi  entanglement entropy [22], Logarithmic Negativity [23], Concurrence [24], Entanglement of formation [25] and Quantum discord [26]. Finally in Sect. 8, we study the concept of non-locality from the violation of Bell-CHSH inequality [27] in De-Sitter space for the two atomic entangled subsystem in the context of OQS.

2 Modelling two atomic open quantum system (OQS)

In nature we always come across open quantum systems where the system interacts with the environment. The most general form of the total Hamiltonian can be written

$$\begin{aligned} H_\mathbf{Total }=H_\mathbf{System }\otimes {\mathcal {I}}_\mathbf{Bath } + {{\mathcal {I}}}_\mathbf{System }\otimes H_\mathbf{Bath }+H_\mathbf{Int }, \end{aligned}$$

(2.1)

where, $H_\mathbf{System }$ represents the two atomic system Hamiltonian,¹ $H_\mathbf{Bath }$ describes the thermal bath Hamiltonian,which is described by massless probe scalar field [28] which is minimally coupled to gravity in static de Sitter background and $H_\mathbf{Int }$ signifies the interaction between the thermal bath and the system under consideration in OQS. ${\mathcal {I}}_\mathbf{System }$ and ${{\mathcal {I}}}_\mathbf{Bath }$ are basically identity operators in system and bath Hilbert space respectively.

In our two atomic OQS set up the system, bath and the interaction Hamiltonian are described by the following expressions [28]:

$$\begin{aligned} H_\mathbf{System }= & {} \frac{\omega }{2}\sum ^{2}_{\alpha =1}{\hat{{\mathbf{n }}}}^{\alpha }.{\sigma }^{\alpha } \end{aligned}$$

(2.2)

$$\begin{aligned} H_\mathbf{Bath }(\tau )= & {} \int ^{\infty }_{0} dr~\int ^{\pi }_{0}d\theta ~\int ^{2\pi }_{0}d\phi ~\left[ \frac{\Pi ^2_{\Phi }(\tau ,r,\theta ,\phi )}{2}\right. \nonumber \\&\quad +\frac{r^2\sin ^2\theta }{2}\left\{ r^2~(\partial _{r}\Phi (\tau ,r,\theta ,\phi ))^2 \right. \nonumber \\&\quad \left. \left. + \frac{\left( (\partial _{\theta }\Phi (\tau ,r,\theta ,\phi ))^2 + \frac{(\partial _{\phi }\Phi (\tau ,r,\theta ,\phi ))^2}{\sin ^2\theta } \right) }{\displaystyle \left( 1-\frac{r^2}{\alpha ^2}\right) }\right\} \right] , \nonumber \\\end{aligned}$$

(2.3)

$$\begin{aligned} H_\mathbf{Int }(\tau )= & {} \mu \sum ^{2}_{\alpha =1}\left( {\hat{\mathbf{n}}}^{\alpha }{\sigma }^{\alpha }\right) \Phi (\tau ,{\mathbf{x}}^{\alpha }) \nonumber \\= & {} \mu \sum ^{2}_{\alpha =1}\left( {\hat{\mathbf{n}}}^{\alpha }.{\sigma }^{\alpha }\right) \Phi (\tau ,{r}^{\alpha },{\theta }^{\alpha },{\phi }^{\alpha }) \end{aligned}$$

(2.4)

It is important to note that, $\omega$ represents the renormalized energy level for two atoms,² given by:

$$\begin{aligned} \begin{array}{lll} \displaystyle \omega =\omega _0+i\times \left\{ \begin{array}{ll} \displaystyle [{\mathcal {K}}^{(11)}(-\omega _{0})-{\mathcal {K}}^{(11)}(\omega _{0})]~~~~~~ &{} {\mathbf{Atom 1 }} \\ \displaystyle [{\mathcal {K}}^{(22)}(-\omega _{0})-{\mathcal {K}}^{(22)}(\omega _{0})]~~~~~~ &{} {\mathbf{Atom 2 }} \end{array} \right. \end{array}\nonumber \\ \end{aligned}$$

(2.5)

Here ${\mathcal {K}}^{\alpha \alpha }(\pm \omega _{0})$ for $\alpha \in \{1,2\}$ are Hilbert transformations of Wightman function computed from the probe massless scalar field, which we have defined explicitly in later section of this paper. Also, $\omega _0$ represents the natural frequency of the two identical atoms, which we have fixed [37,38,39,40,41] at:³

$$\begin{aligned}&\omega _0=\frac{i}{k}\left( n+\frac{1}{2}\right) ~~~\forall ~n\in {\mathbb {Z}}~~~\mathrm{and} \nonumber \\&k=\sqrt{\alpha ^2-r^2}>0, ~\alpha =\sqrt{\frac{3}{\Lambda }}>0~~\mathrm{as}~~\Lambda >0 \end{aligned}$$

(2.6)

for rest of the computation performed in this paper. In this context, the atoms are characterised by the label $\alpha \in [1,2]$ and $\sigma ^{\alpha }_{i}\forall i \in [1,2,3]$ are the Pauli spin matrices.

The bath Hamiltonian for the probe massless scalar field have been expressed in the static patch of de Sitter space. The prime reason for choosing the de Sitter background can be understood from the fact that it is the simplest non trivial curved spacetime and is maximally symmetric. This unique feature of the de Sitter spacetime has led to extensive study on the quantization of the scalar fields [29,30,31,32,33,34] in this spacetime, which actually acts as the bath for our model. Apart from this the significance of this spacetime can also be understood in the context of cosmology. Current cosmological observations in association with the inflation theory suggets that our universe approches the exponentially expanding de Sitter phase both in the early and the late time scales ,i.e in both the far past and the far future.

In static patch of de Sitter space the background space time metric is described by the following infinitesimal line element:

$$\begin{aligned} ds^{2}= & {} \left( 1-\frac{r^{2}}{\alpha ^{2}}\right) dt^{2} -\left( 1-\frac{r^{2}}{\alpha ^{2}}\right) ^{-1}dr^{2} \nonumber \\&-r^{2}(d\theta ^{2}+\sin ^{2}\theta d\phi ^{2})~~~\mathrm{where}~~\alpha =\sqrt{\frac{3}{\Lambda }}>0 \end{aligned}$$

(2.8)

Instead of using the time variable as t in the present context we have introduced a new rescaled time variable $\tau$, which is defined as:

$$\begin{aligned}&\tau =\sqrt{g_{00}}~t=\frac{k}{\alpha }~t=\sqrt{1-\frac{r^2}{\alpha ^2}}~t \nonumber \\&\quad \mathrm{where}~~k=\sqrt{\alpha ^2-r^2}>0 \end{aligned}$$

(2.9)

Under this assumption of background space time, the massless probe scalar field is minimally coupled to the gravity in this context. Additionally, it is important to mention here that the interaction is controlled by the interaction strength or coupling parameter $\mu$ through which the bath degrees of freedom is coupled to the two atomic subsystem. For the description of the OQS we will use the traditional approach of solving the GSKL Master Equation which describes the non-unitary time evolution of the reduced density matrix of the subsystem.

3 Non unitary time evolution of the reduced subsystem

The prime objective in this section is to describe the time evolution of an OQS and to arrive at GSKL master equation which properly describes non-unitary behaviour and can be obtained by performing partial trace over the bath content i.e. the massless probe scalar field placed at the static patch of the de Sitter background space time.

The exact mathematical form of the quantum correlation part of the total density matrix can only be obtained by solving the GSKL master equation. The overall time evolution of the total density matrix is given by the following Liouville Von Neumann equation.

$$\begin{aligned} \frac{d}{d\tau }\rho _\mathbf{Total }(\tau )=-i [H_\mathbf{Total }, \rho _\mathbf{Total }(\tau )] \end{aligned}$$

(3.1)

One can perform a unitary transformation to bring the above Liouville Von Neumann equation into the following convenient form

$$\begin{aligned} \frac{d}{d\tau }\rho _\mathbf{Total }(\tau )=-i [H_\mathbf{Int }(\tau ), \rho _\mathbf{Total }(\tau )] \end{aligned}$$

(3.2)

We can integrate the above equation to obtain

$$\begin{aligned} \rho _\mathbf{Total }(\tau ) = \rho _\mathbf{Total }(0) - i \int \limits _0^t ds [H_\mathbf{Int }(s),\rho _\mathbf{Total }(s)] \end{aligned}$$

By inserting this back into Eq. (3.2) and tracing out the bath degrees of freedom, we arrive at

$$\begin{aligned}&\frac{d}{dt}\rho _\mathbf{System }(\tau ) \nonumber \\&\quad = -\int \limits _0^t ds {\text {Tr}}_\mathbf{Bath } \{[H_\mathbf{Int }(t), [H_I(s), \rho _\mathbf{Total }(s)]] \} \end{aligned}$$

(3.3)

where we have taken

$$ {\text {Tr}}_B\{[H_\mathbf{Int }(t),\rho _\mathbf{Total }(0)]\} = 0 $$

which means that initially the interaction does not create any dynamics in the bath.

Equation (3.3) is the integro-differential equation which is non-Markovian in nature as the dynamics at a particular time depends on its past evolution. This is significantly difficult to calculate. Besides on the right hand side of Eq. (3.3) we still have the total density operator.

Since, we are primarily concerned about the two atomic system, by doing the following approximations we bring this non-Markovian integro-differential equation into a Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) [19, 20] master equation which is Markovian in nature and can be solved easily.

• Born Approximation: Generally in the treatment of any open quantum system the density matrix is usually represented by

  $$ \rho _\mathbf{Total }(t)=\rho _\mathbf{System }(t)\otimes \rho _\mathbf{Bath }(t)+\rho _\mathbf{Correlation }(t) $$

  We assume that the interaction between the bath and the system is switched on at time t=0, which is a good enough assumption because finding such a time scale is always posiible before which there was no interaction i.e ($\rho _\mathbf{Correlation }(0)=0$) . Besides, we make a stronger claim that the coupling between the system and the bath is weak such that the influence of the bath is small. This is known as Born approximation. This approximation is important since we can ignore the initial corelation between the atomic system and bath. In this weak coupling limit situation, the time evolved density matrix is reduced to the form

  $$ \rho _\mathbf{Total }(t) \approx \rho _\mathbf{System }(t) \otimes \rho _\mathbf{Bath }(t) $$

  This assumption reduces the Master equation to Redfield master equation which is still non-Markovian in nature.

• Markov Approximation: To obtain a Markovian master equation, we assume the relaxation time($\tau _{R}$) to be more compared to the bath correlation time($\tau _{B}$), i.e, $\displaystyle \tau _{R}\gg \tau _{B}$

  With this approximation we may take the bath as nearly constant i.e. $\displaystyle \rho _\mathbf{Bath }(t) \approx \rho _\mathbf{Bath }(0)$

  As a result the total density operator can be simplified to

  $$ \rho _\mathbf{Total }(t) = \rho _\mathbf{System }(t) \otimes \rho _\mathbf{Bath }(0) $$

  This approximation grouped with the Born approximation is often regarded as the Born-Markov approximation. However, with this approximation alone the resulting master equation does not guarantee to generate a quantum synamical semigroup. We need one more approximation to write the Master equation in Lindblad form.

• Secular Approximation: With this approximation we can average out and discard highly oscillating terms, comapred to the system timescale of interest, in the Markovian master equation.

With the aid of all these approximations we obtain the following markovian master equation

$$\begin{aligned}&\frac{d}{d\tau }\rho _\mathbf{System }(\tau ) \nonumber \\&\quad =-i[H_\mathbf{eff }(\tau ),\rho _\mathbf{System }(\tau )]+{\mathcal {L}} [\rho _\mathbf{System }(\tau )] \end{aligned}$$

(3.4)

where $H_\mathbf{eff }$ is the effective Hamiltonian and ${\mathcal {L}} [\rho _\mathbf{System }(\tau )]$ is the quantum dissipator or the Lindbladian operator. These are discussed in more detail in the following sections.

4 Effective Hamiltonian construction

For our system, the effective Hamiltonian can be expressed as

$$\begin{aligned} H_\mathbf{eff }= & {} H_\mathbf{System }+H_\mathbf{Lamb~ Shift } \nonumber \\= & {} \underbrace{\frac{\omega }{2}\sum _{\alpha =1}^2 {\mathbf{n}}^{\alpha }.{\sigma }^{\alpha }}_{\mathbf{Two Atomic System}}-\underbrace{\frac{i}{2}\sum _{\alpha \beta =1}^{2} \sum _{ij=1}^3 H_{ij}^{\alpha \beta }({\mathbf{n }_i}^{\alpha }.{\sigma _i}^{\alpha })({\mathbf{n }_j}^{\alpha }.{\sigma _j}^{\alpha })}_{\mathbf{Lamb Shift=Heisenberg spin chain}},\nonumber \\ \end{aligned}$$

(4.1)

where the first term in the effective Hamiltonian, physically represents the Hamiltonian of the two atomic system whereas the second term is known as the Lamb Shift Hamiltonian [16, 46], which characterizes the atomic Lamb Shift that occurs due to the interaction between the massless free probe scalar field with the two atomic system under consideration in the background of static de Sitter space. It actually measures a shift in the energy levels due to this interaction.

Here ${\mathbf{n}^\alpha }$ and ${\mathbf{n}^\beta }$ represent the normal unit vectors of the two atoms under consideration in the present OQS set up. The angles between the normal unit vectors and Pauli spin matrices are characterized by the three Euler Angles $\alpha$,$\beta$ and $\gamma$. However, we consider here that these Euler angles for two atoms are different to get more general result.

Therefore the Lamb Shift Hamiltonian can be re-expressed in terms of Euler angles as:

$$\begin{aligned} H_\mathbf{Lamb~ Shift }= & {} -\frac{i}{2}\sum _{\alpha \beta =1}^{2} \sum _{ij=1}^3 H_{ij}^{\alpha \beta }({\mathbf{n }_i}^{\alpha }.{\sigma _i}^{\alpha })({\mathbf{n }_j}^{\alpha }.{\sigma _j}^{\alpha }) \nonumber \\= & {} -\frac{i}{2}\sum _{\alpha \beta =1}^{2} \sum _{ij=1}^3 H_{ij}^{\alpha \beta } \cos (\alpha _i^{\alpha }) \cos (\alpha _j^{\beta })\sigma _i^{\alpha }\sigma _j^{\beta }.\nonumber \\ \end{aligned}$$

(4.2)

In this present context, we define the following sets of Pauli operators for the two atoms using tensor product:

$$\begin{aligned} \sigma _i^1= & {} \sigma _i \otimes \sigma _0 ~~~~~(\mathbf{Atom 1}) ,~~~~~~~ \sigma _i^2=\sigma _0 \otimes \sigma _i ~~~~~(\mathbf{Atom 2})\nonumber \\ \end{aligned}$$

(4.3)

Here $\sigma _0$ is the $2 \times 2$ identity matrix and $\sigma _i$ ($i=1,2,3$) are usual Pauli matrices. In our case, we have changed the basis for representing the effective Hamiltonian from the $\sigma _1$,$\sigma _2$,$\sigma _3$ to the $\sigma _+$,$\sigma _-$,$\sigma _3$ basis which along with the identity matrix forms a complete basis for the space of $2 \times 2$ matrices. The change of basis basically reduces the number of differential equations that is obtained from the GSKL master equation and makes them simpler to solve.

In the transformed basis $\sigma _+$ and $\sigma _-$ are defined as: follows:

$$\begin{aligned} \sigma _+= & {} \frac{1}{2}(\sigma _1+i \sigma _2) =\begin{pmatrix} 0 &{} &{} &{} &{} &{} 1\\ 0 &{} &{} &{} &{} &{}0 \end{pmatrix}, \nonumber \\ \sigma _-= & {} \frac{1}{2}(\sigma _1-i \sigma _2) =\begin{pmatrix} 0 &{} &{} &{} &{} &{} 0\\ 1 &{} &{} &{} &{} &{} 0 \end{pmatrix}. \end{aligned}$$

(4.4)

Therefore the operators $\sigma _+^1$,and $\sigma _-^1$ for Atom 1 is defined as:

$$\begin{aligned} \sigma _+^1= & {} (\sigma _+ \otimes \sigma _0) =\begin{pmatrix} 0 &{} &{} &{} &{} &{} \sigma _0\\ 0 &{} &{} &{} &{} &{} 0 \end{pmatrix}, \nonumber \\ \sigma _-^1= & {} (\sigma _-\otimes \sigma _0) =\begin{pmatrix} 0 &{} &{} &{} &{} &{} 0\\ \sigma _0 &{} &{} &{} &{} &{} 0 \end{pmatrix} \end{aligned}$$

(4.5)

Similarly, the operators $\sigma _+^2$,and $\sigma _-^2$ for Atom 2 is defined in a similar way as:

$$\begin{aligned} \sigma _+^2= & {} (\sigma _0 \otimes \sigma _+) =\begin{pmatrix} \sigma _+ &{} 0\\ 0 &{} \sigma _+ \end{pmatrix}, \nonumber \\ \sigma _-^2= & {} (\sigma _0\otimes \sigma _-) =\begin{pmatrix} \sigma _- &{} &{} 0\\ 0 &{} &{} \sigma _- \end{pmatrix} \end{aligned}$$

(4.6)

Now we transform the Hamiltonian in the new basis which is equivalent to diagonalising the $n_i^{\alpha }\sigma _i^{\beta }$ term, and basically reduces it to the $\sigma _3$ model. Hence the terms of $H_\mathbf{System }$ in the diagonalized basis respectively becomes:

$$\begin{aligned} \mathbf{Atom~1:}~~~ H_1= & {} \frac{\omega }{2} n_i^1\sigma _i^1=\frac{\omega }{2} n_i\left( \sigma _i \otimes \sigma _0\right) \Longrightarrow \nonumber \\&\frac{\omega '}{2} \sigma _3 \otimes \sigma _0 \end{aligned}$$

(4.7)

$$\begin{aligned} \mathbf{Atom~2:}~~~ H_2= & {} \frac{\omega }{2} n_i^2\sigma _i^2=\frac{\omega }{2} n_i\left( \sigma _0 \otimes \sigma _i \right) \Longrightarrow \nonumber \\&\frac{\omega '}{2} \sigma _0 \otimes \sigma _3 \end{aligned}$$

(4.8)

where $\omega '$ is the modified frequency. When the term $n_i\sigma _i$ is diagonalized in the new basis, a factor of $\sqrt{n_3^2+n_+n_{-}}$ arises, which can be incorporated in the frequency as the modification factor. Thus the frequency gets modified by this process and is given by:

$$\begin{aligned} \omega '=\omega \sqrt{n_3^2+n_+n_-}. \end{aligned}$$

(4.9)

Here, $n_3$, $n_+$ and $n_-$ are the normal unit vectors of the two atoms in the new basis defined as:

$$\begin{aligned} n_+= & {} \frac{1}{2}(n_1+i n_2)=\frac{1}{2}(\cos \alpha _1+ i\cos \alpha _2) , \nonumber \\ n_-= & {} \frac{1}{2}(n_1-i n_2)=\frac{1}{2}(\cos \alpha _1-i\cos \alpha _2) \end{aligned}$$

(4.10)

In the new diagonalized basis (already mentioned earlier) the above mentioned Lamb Shift Hamiltonian part reduces to the following simplified form:

$$\begin{aligned} H_\mathbf{Lamb Shift }=-\frac{i}{2}\sum _{\alpha , \beta =1}^2 \sum _{i,j=\pm }^3 H_{ij}^{\alpha \beta } \sigma _i^\alpha \sigma _j^\beta \end{aligned}$$

(4.11)

In this context, the effective Hamiltonian matrix elements, $H_{ij}^{\alpha \beta }$, is given by:

$$\begin{aligned} H_{ij}^{\alpha \beta }={ A^{\alpha \beta }}\delta _{ij}-i { B^{\alpha \beta }}\epsilon _{ijk}\delta _{3j}-{{ A^{\alpha \beta }}\delta _{3i}\delta _{3j}} \end{aligned}$$

(4.12)

where $A^{\alpha \beta }$ and $B^{\alpha \beta }$ for the two atomic system are defined as :

$$\begin{aligned} A^{\alpha \beta }= & {} \frac{\mu ^2}{4}[K^{\alpha \beta }(\omega _0)+K^{\alpha \beta }(-\omega _0)], \nonumber \\ B^{\alpha \beta }= & {} \frac{\mu ^2}{4}[K^{\alpha \beta }(\omega _0)-K^{\alpha \beta }(-\omega _0)], \end{aligned}$$

(4.13)

where $K^{\alpha \beta }$ is basically the Hilbert Transform of the Wightman function (two point correlator) computed from the probe massless scalar field placed at the bath and is given by the following expression:

$$\begin{aligned} K^{\alpha \beta }(\pm \omega _0)= & {} \frac{P}{\pi i}\int _{-\infty }^{\infty }d\omega \frac{{\mathcal {G}}^{\alpha \beta }(\pm \omega )}{\omega \pm \omega _0} \nonumber \\= & {} \frac{P}{\pi i}\int _{-\infty }^{\infty }d\omega \frac{1}{\omega \pm \omega _0} \int _{-\infty }^{\infty }d\Delta \tau e^{\pm i \omega \Delta \tau } G^{\alpha \beta }(\Delta \tau )\nonumber \\ \end{aligned}$$

(4.14)

where P is the principal value of the integral. Here ${\mathcal {G}}^{\alpha \beta }$ is the Fourier transform of the of the Wightman function in the frequency ($\omega$) space and can be expressed as:

$$\begin{aligned} {\mathcal {G}}^{\alpha \beta } (\pm \omega _0)=\int _{-\infty }^{\infty }d\Delta \tau ~ e^{\pm \iota \omega \Delta \tau }~G^{\alpha \beta }(\Delta \tau ), \end{aligned}$$

(4.15)

where $\omega _0$ is the energy difference between the ground and excited states of the atoms and $G^{\alpha \beta }$ is the forward two atomic Wightman Function which is defined as:⁴

$$\begin{aligned}&G^{\alpha \gamma }(\Delta \tau =\tau -\tau ^{\prime }) \nonumber \\&\quad =\langle \Phi (x_\alpha ,\tau ) \Phi (x_\gamma ,\tau ^{\prime })\rangle _{\beta } \nonumber \\&\quad = \frac{\int _\mathbf{S {} \mathbf{K} } {{\mathcal {D}}}\Psi ~\langle \Psi |\exp (-\beta ~H_\mathbf{bath })~\Phi (x_\alpha ,\tau ) \Phi (x_\gamma ,\tau ^{\prime })|\Psi \rangle }{\int _\mathbf{SK } {\mathcal {D}} \Psi ~\langle \Psi |\exp (-\beta ~H_\mathbf{bath })|\Psi \rangle }.\nonumber \\ \end{aligned}$$

(4.16)

In the above equations $\mu$ the coupling parameter, represents the interaction strength between the system and the external thermal bath (i.e the gravitationally coupled scalar) field degrees of freedom. The structure of the elements of the coefficient matrix $H_{ij}^{\alpha \beta }$ can be computed in terms of the Wightman function of the external free probe massless scalar field in static De-Sitter background, which finally fix the structure of the effective Hamiltonian in the present scenario.

5 Quantum dissipator or Lindbladian construction

The concept of fluctuation and dissipation in the context of OQS is introduced into the system by the additional contribution of the Lindbladian operator in the time evolution equation of the reduced subsytem density matrix. The second term in the Gorini Kossakowski Sudarshan Lindblad Master (GSKL) equation is actually characterised as the Lindbladian or Quantum Dissipator which in our present context can be written as:

$$\begin{aligned} \mathcal {L}[\rho _\mathbf{System }(\tau )]= & {} \frac{1}{2}\sum _{i,j=1}^3 \sum _{\alpha ,\beta =1}^2 C_{ij}^{\alpha \beta } \nonumber \\&\left[ 2(n_j^\beta \cdot \sigma _j^\beta )\rho _\mathbf{System (\tau )}(n_i^\alpha \cdot \sigma _i^\alpha ) \right. \nonumber \\&\qquad \left. -\left\{ (n_i^\alpha \cdot \sigma _i^\alpha )(n_j^\beta \cdot \sigma _j^\beta ),\rho _\mathbf{System }(\tau )\right\} \right] ,\nonumber \\ \end{aligned}$$

(5.1)

where $\rho _\mathbf{System }$ is the reduced subsystem density matrix of the two entangled atomic system obtained after partially tracing over the external bath scalar field degrees of freedom. The coefficient matrix $C_{ij}^{\alpha \beta }$ is known as the Gorini Kossakowski Sudarshan Lindblad (GSKL) matrix, which is constructed under the weak coupling limiting approximation on the coupling parameter $\mu$ as appeared in the interaction Hamiltonian. In the context of OQS, the Lindbladian captures the effect of dissipation.

In the transformed basis, the Lindbladian can be re-expressed as:

$$\begin{aligned} \mathcal {L}[\rho _\mathbf{System }(\tau )]= & {} \frac{1}{2}\sum _{i,j=\pm }^3 \sum _{\alpha , \beta =1}^2 C_{ij}^{\alpha \beta } \nonumber \\&\left[ 2 \sigma _j^{\beta }\rho _\mathbf{System }(\tau ) \sigma _i^\alpha \right. \nonumber \\&\quad \left. -\left\{ \sigma _i^\alpha \sigma _j^\beta ,\rho _\mathbf{System }(\tau )\right\} \right] . \end{aligned}$$

(5.2)

The matrix GSKL matrix $C_{ij}^{\alpha \beta }$ is given by the following expression:

$$\begin{aligned} C_{ij}^{\alpha \beta }={\tilde{A}}^{\alpha \beta } \delta _{ij}-i {\tilde{B}}^{\alpha \beta }\epsilon _{ijk}\delta _{3k}-{\tilde{A}}^{\alpha \beta }\delta _{3k}\delta _{3j} \end{aligned}$$

(5.3)

where the quantities ${\tilde{A}}^{\alpha \beta }$ and ${\tilde{B}}^{\alpha \beta }$ for the two atomic system is defined as:

$$\begin{aligned} {\tilde{A}}^{\alpha \beta }= & {} \frac{\mu ^2}{4}\left[ {\mathcal {G}}^{\alpha \beta } (\omega _0)+{\mathcal {G}}^{\alpha \beta } (-\omega _0)\right] \nonumber \\= & {} \frac{\mu ^2}{4}\int _{-\infty }^{\infty }G^{\alpha \beta }(\Delta \tau )\left[ e^{\iota \omega _0 \Delta \tau }+e^{-\iota \omega _0 \Delta \tau }\right] \end{aligned}$$

(5.4)

$$\begin{aligned} {\tilde{B}}^{\alpha \beta }= & {} \frac{\mu ^2}{4}\left[ {\mathcal {G}}^{\alpha \beta }(\omega _0)-{\mathcal {G}}^{\alpha \beta }(-\omega _0)\right] \nonumber \\= & {} \frac{\mu ^2}{4}\int _{-\infty }^{\infty }G^{\alpha \beta }(\Delta \tau )\left[ e^{i \omega _0 \Delta \tau }-e^{-i \omega _0 \Delta \tau }\right] \end{aligned}$$

(5.5)

The components of $C_{ij}^{\alpha \beta }$ matrix are given in the Appendix B.

6 Time evolution of the reduced subsystem density matrix

The density matrices for Atom 1 and Atom 2 are given by following expressions in the Bloch sphere representation:

$$\begin{aligned}&\mathbf{Atom 1:} \nonumber \\&\rho _1(\tau ) =\frac{1}{2}\left( I+\sum _{i=1}^3 a_i (\tau )\sigma _i \right) =\frac{1}{2}\left[ I+\mathbf{a}(\tau ) \cdot {\sigma } \right] . \end{aligned}$$

(6.1)

$$\begin{aligned}&\mathbf{Atom 2:} \nonumber \\&\rho _2(\tau )=\frac{1}{2}\left( I+\sum _{j=1}^3 b_j(\tau ) \sigma _j \right) =\frac{1}{2}\left[ I+\mathbf{b}(\tau ) \cdot {\sigma } \right] . \end{aligned}$$

(6.2)

Since the two atoms are initially not entangled, the density matrix for the system can be written as the direct product of the two individual density matrices i.e.

$$\begin{aligned}&\rho _\mathbf{System }(\tau ) \nonumber \\&\quad =\rho _1(\tau )\otimes \rho _2(\tau )=\frac{1}{2} \left( I+\sum _{i=1}^3 a_i(\tau ) \sigma _i \right) \nonumber \\&\qquad \otimes \frac{1}{2}\left( I+\sum _{j=1}^3 b_j(\tau ) \sigma _j \right) \nonumber \\&\quad = \frac{1}{4}\left[ \sigma _0\otimes \sigma _0 +\sum _{j=1}^3 b_j(\tau ) \sigma _0 \otimes \sigma _j \right. \nonumber \\&\qquad \left. + \sum _{i=1}^3 a_i(\tau ) \sigma _i \otimes \sigma _0+\sum _{i,j=1}^3 a_i(\tau ) b_j(\tau ) \sigma _i \otimes \sigma _j \right] .\nonumber \\ \end{aligned}$$

(6.3)

Note that in the above equation the identity matrix is denoted by $\sigma _0$. Now we define:

$$\begin{aligned}&a_i(\tau )\equiv a_{i0}(\tau )~~~~~~~~~~~\forall ~~i=1,2,3, \end{aligned}$$

(6.4)

$$\begin{aligned}&b_j(\tau )\equiv a_{0j}(\tau )~~~~~~~~~~~\forall ~~j=1,2,3, \end{aligned}$$

(6.5)

$$\begin{aligned}&a_i(\tau ) b_j(\tau )\equiv a_{ij}(\tau )~~~~~\forall ~~i,j=1,2,3. \end{aligned}$$

(6.6)

Using these definitions the density matrix $\rho _\mathbf{System }(\tau )$ for the reduced subsystem can be re-expressed as:

$$\begin{aligned}&\rho _\mathbf{System }(\tau ) \nonumber \\&\quad =\frac{1}{4}\left[ \sigma _0 \otimes \sigma _0 +\sum _{i=1}^3 a_{0i}(\tau )(\sigma _0 \otimes \sigma _i) \right. \nonumber \\&\qquad \left. +\sum _{i=1}^3 a_{i0}(\tau )(\sigma _i \otimes \sigma _0)+\sum _{i,j=1}^3 a_{0i}(\tau )(\sigma _i \otimes \sigma _j)\right] \end{aligned}$$

(6.7)

In the new transformed basis i.e in terms of $\sigma _+$, $\sigma _-$ and $\sigma _3$ basis $\rho _\mathbf{System }(\tau )$ can be written as:

$$\begin{aligned}&\rho _\mathbf{System }(\tau ) \nonumber \\&\quad =\frac{1}{4}\left[ \sigma _0 \otimes \sigma _0 +\sum _{m=+,-}^3 a_{0m}(\tau )(\sigma _0 \otimes \sigma _m) \right. \nonumber \\&\qquad \left. +\sum _{m=+,-}^3 a_{m0}(\tau )(\sigma _m \otimes \sigma _0)+\sum _{m,n=+,-}^3 a_{mn}(\tau )(\sigma _i \otimes \sigma _j)\right] .\nonumber \\ \end{aligned}$$

(6.8)

In this new basis the reduced subsystem density matrix, in terms of the Bloch vectors can be explicitly written as:

$$\begin{aligned}&\rho _\mathbf{System }(\tau )=\frac{1}{4}\begin{pmatrix} 1+a_{03}+a_{30}+a_{33} &{} 0 &{} 0 &{} a_{++}\\ 0 &{} 1-a_{03}+a_{30}-a_{33} &{} a_{+-} &{} 0\\ 0 &{} a_{-+} &{} 1+a_{03}-a_{30}-a_{33} &{} 0\\ a_{--} &{} 0 &{} 0 &{} 1-a_{03}-a_{30}+a_{33} \end{pmatrix}, \end{aligned}$$

(6.9)

where the Hermiticity and $\mathrm{Tr}(\rho _\mathbf{System }(\tau ) )=1$ property of the density matrix has been used to find the matrix elements explicitly.

6.1 Large scale time dependent solution

Solving the master equation in the new basis i.e. $\sigma _+$, $\sigma _-$ and $\sigma _3$ basis in the large time limit ($\tau =\infty$) we get the following solution for the components of the density matrix as given by:

$$\begin{aligned} \begin{aligned} a_{03}(\infty )=a_{30}(\infty )= & {} -\mathrm{tanh(\pi k \omega )} \\ a_{33}(\infty )= & {} \mathrm{tanh^2(\pi k \omega )} \\ a_{++}(\infty )=a_{--}(\infty )= & {} 0 \\ a_{+-}(\infty )=a_{-+}(\infty )= & {} 0 \end{aligned} \end{aligned}$$

(6.10)

where all other Bloch vector components are zero.

We can get the large time reduced density matrix from the above solution, which can be expressed as

$$\begin{aligned}&\rho _\mathbf{System }(\infty ) \nonumber \\&\quad = \frac{1}{4}\left[ \sigma _0 \otimes \sigma _0 + a_{03}(\infty )(\sigma _0 \otimes \sigma _3) \right. \nonumber \\&\qquad \left. +a_{30}(\infty )(\sigma _3 \otimes \sigma _0)+ a_{33}(\infty )(\sigma _3 \otimes \sigma _3)\right] \nonumber \\&\quad = \frac{1}{4}\begin{pmatrix} 1+2a_{03}(\infty )+a_{33}(\infty ) &{} &{} 0 &{} &{} 0 &{} &{} 0\\ 0 &{} &{} 1-a_{33}(\infty ) &{} &{} 0 &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} 1-a_{33}(\infty ) &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} 0 &{} &{} 1-2a_{03}(\infty )+a_{33}(\infty ) \end{pmatrix} \end{aligned}$$

(6.11)

The large time behaviour basically demonstrates the equilibrium behaviour of the system. Hence the solution of the Bloch vectors obtained in the large time scale is applicable in any basis. These solutions can therefore be used as the boundary conditions for obtaining the finite time solution of the density matrix. Writing the density matrix in terms of the solutions of the Bloch vectors as,

$$\begin{aligned}&\rho _\mathbf{System }(\infty ) =\frac{1}{4}\begin{pmatrix} (1-\mathrm{tanh(\pi k \omega )})^2 &{} &{} 0 &{} &{} 0 &{} &{} 0\\ 0 &{} &{} 1-\mathrm{tanh^2(\pi k \omega )} &{} &{} 0 &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} \mathrm{1-tanh^2(\pi k \omega )} &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} 0 &{} &{} \mathrm{(1+tanh(\pi k \omega ))^2} \end{pmatrix} \end{aligned}$$

(6.12)

On the other hand, from quantum statistical mechanics one can compute the expression for the density matrix at finite temperature and large time limit, which is given by the following expression:

$$\begin{aligned}&\rho _\mathbf{System }(\infty ) = \frac{e^{-\beta H_\mathbf{System }}}{\mathrm{Tr} (e^{-\beta H_\mathbf{System }})}\nonumber \\&\quad =\frac{1}{4}\begin{pmatrix} \left( 1-\mathrm{tanh\left( \frac{\beta \omega }{2}\right) }\right) ^2 &{} &{} 0 &{} &{} 0 &{} &{} 0\\ 0 &{} &{} 1-\mathrm{tanh^2\left( \frac{\beta \omega }{2}\right) } &{} &{} 0 &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} \mathrm{1-tanh^2\left( \frac{\beta \omega }{2}\right) } &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} 0 &{} &{} \mathrm{\left( 1+tanh\left( \frac{\beta \omega }{2}\right) \right) ^2} \end{pmatrix}. \end{aligned}$$

(6.13)

Comparing Eq (6.12) and Eq (6.13), we obtain the following result:

$$\begin{aligned}&T=\frac{1}{\beta }=\frac{1}{2\pi k}=\frac{1}{2 \pi \sqrt{\alpha ^2-r^2}}, ~~~\mathrm{where} \nonumber \\&\quad \alpha =\sqrt{\frac{3}{\Lambda }}>0. \end{aligned}$$

(6.14)

which physically represents the equilibrium temperature of the thermal bath at large time scale.

In ref. [28] we have explicitly shown that the temperature of the thermal bath can be computed in terms of the Gibbons Hawking temperature and Unruh temperature, which can be expressed as [49]:

$$\begin{aligned}&T=\sqrt{T^2_\mathbf{GH }+T^2_\mathbf{Unruh }}=\frac{1}{2\pi \alpha }\sqrt{1+\frac{r^2}{\alpha ^2-r^2}}, ~~~\mathrm{where} \nonumber \\&\quad \alpha =\sqrt{\frac{3}{\Lambda }}>0. \end{aligned}$$

(6.15)

Here Gibbons Hawking and Unruh temperature is defined in the present context as:

$$\begin{aligned}&T_\mathbf{GH }=\frac{1}{2\pi \alpha },~~ T_\mathbf{Unruh }=\frac{a}{2\pi }=T_\mathbf{GH }\frac{r}{\sqrt{\alpha ^2-r^2}}, ~~\mathrm{where} \nonumber \\&\quad a=2\pi T_\mathbf{GH }\frac{r}{\sqrt{\alpha ^2-r^2}}=\frac{1}{\alpha }\frac{r}{\sqrt{\alpha ^2-r^2}}.\end{aligned}$$

(6.16)

Now, we know that in static patch of the de Sitter space [35, 36] the curvature is determined by the following expression for the Ricci scalar:

$$\begin{aligned} R=\frac{12}{\alpha ^2}>0,~~~\mathrm{where}~~\alpha =\sqrt{\frac{3}{\Lambda }}>0. \end{aligned}$$

(6.17)

Consequently, the equilibrium temperature of the thermal bath can be re-expressed in terms of the curvature of the static patch of the de Sitter space as:

$$\begin{aligned} T=\frac{1}{\beta }=\frac{1}{2\pi k}=\frac{1}{2 \pi \sqrt{\left( \frac{12}{R}\right) -r^2}}=\frac{\sqrt{R}}{2\pi }\frac{1}{\sqrt{12-Rr^2}}.\nonumber \\ \end{aligned}$$

(6.18)

Similarly, the Gibbons Hawking and Unruh temperature can be re-expressed in terms of the curvature of the static patch of the de Sitter space as:

$$\begin{aligned} T_\mathbf{GH }= & {} \frac{\sqrt{R}}{4\sqrt{3}\pi },~~~~~ T_\mathbf{Unruh }=\frac{a}{2\pi }=T_\mathbf{GH }\frac{\sqrt{R}r}{\sqrt{12-Rr^2}}.\nonumber \\ \end{aligned}$$

(6.19)

In this context, one can consider the following limiting situations:

• 1. Flat space limit:

  Flat space limit is characterised by the following condition:

  $$\begin{aligned}&R\rightarrow 0,\Longrightarrow \alpha \rightarrow \infty ,\Longrightarrow T_\mathbf{GH }, \nonumber \\&T_\mathbf{Unruh }\rightarrow 0,\Longrightarrow T\rightarrow 0,\Longrightarrow k\rightarrow \infty \Longrightarrow k>>L.\nonumber \\ \end{aligned}$$

  (6.20)

  Here L represents the Euclidean distance between the two atoms. In this case, we will get back the result obtained in the Minkowski flat space inertial case where $k>>L$.

• 2. Zero acceleration limit:

  The zero acceleration limit is characterised by the following condition:

  $$\begin{aligned}&a\rightarrow 0,\Longrightarrow r\rightarrow 0,\Longrightarrow T_\mathbf{GH }\ne 0, \nonumber \\&T_\mathbf{Unruh }\rightarrow 0,\Longrightarrow T\rightarrow T_\mathbf{GH },\Longrightarrow k<<L. \end{aligned}$$

  (6.21)

  In this case, we will get back the result obtained in the limiting case where $k<<L$.

6.2 Arbitrary time dependent general solution

To obtain the finite time solution of the Bloch vector components, we use the $\sigma _+$, $\sigma _-$ and $\sigma _3$ basis. Substituting the components of the density matrix in the GSKL master equation, in this new basis, we obtain the evolution equations for the Bloch vector components.

$$\begin{aligned} {\dot{a}}_{03}(\tau )= & {} \frac{1}{4}(A^{12}+A^{21})(a_{++}-a_{--}) \nonumber \\&+\frac{1}{2}({\tilde{A}}^{21}-{\tilde{A}}^{12})(a_{++}-a_{--})\nonumber \\&+\frac{i}{2}({\tilde{B}}^{12}+{\tilde{B}}^21)(a_{+-}+a_{-+})+4i {\tilde{B}}^{22}\nonumber \\ {\dot{a}}_{03}(\tau )= & {} \frac{1}{4}(A^{12}+A^{21})(a_{++}-a_{--})\nonumber \\&+\frac{1}{2}({\tilde{A}}^{21}+{\tilde{A}}^{12})(a_{++}-a_{--})\nonumber \\&+\frac{i}{2}({\tilde{B}}^{12}+{\tilde{B}}^{21})(a_{+-}+a_{-+}) +4i {\tilde{B}}^{11}\nonumber \\ {\dot{a}}_{++}(\tau )= & {} (A^{12}+A^{21})(a_{03}+a_{30})+i a_{++}(B^{11}+B^{22})\nonumber \\&+2\omega a_{++}+2 {\tilde{A}}^{22}a_{+-}+2 {\tilde{A}}^{11}a_{-+}\nonumber \\&+2{\tilde{A}}^{21}(a_{03}-a_{30}-2a_{33})+2{\tilde{A}}^{12}(-a_{03}+a_{30}-2a_{33}) \nonumber \\ {\dot{a}}_{+-}(\tau )= & {} i (-B^{12}+B^{21})(a_{30}-a_{03})+i a_{12}(B^{11}-B^{22})\nonumber \\&+2i {\tilde{B}}^{21}(a_{03}-a_{30}+2a_{33})\nonumber \\&-2i {\tilde{B}}^{12}(a_{03}-a_{30}+2a_{33})+2{\tilde{A}}^{11}a_{--}+2{\tilde{A}}^{22}a_{++} \nonumber \\ {\dot{a}}_{-+}(\tau )= & {} i (B^{12}-B^{21})(a_{03}-a_{30})+i a_{21}(-B^{11}+B^{22})\nonumber \\&-2i {\tilde{B}}^{21}(a_{03}+a_{30}+2a_{33})\nonumber \\&-2i {\tilde{B}}^{12}(-a_{03}-a_{30}+2a_{33})+2{\tilde{A}}^{11}a_{++}+2{\tilde{A}}^{22}a_{--} \nonumber \\ {\dot{a}}_{--}(\tau )= & {} (A^{12}+A^{21})(-a_{03}-a_{30})+i a_{--}(B^{11}+B^{22})\nonumber \\&+2\omega a_{--}+2 {\tilde{A}}^{22}a_{-+}+2 {\tilde{A}}^{11}a_{+-}\nonumber \\&+2{\tilde{A}}^{21}(a_{03}-a_{30}-2a_{33})+2{\tilde{A}}^{12}(-a_{03}+a_{30}-2a_{33}) \nonumber \\ {\dot{a}}_{33}(\tau )= & {} -({\tilde{A}}^{12}+{\tilde{A}}^{21})a_{++}\nonumber \\&-({\tilde{A}}^{21}+{\tilde{A}}^{12})a_{--}+4i {\tilde{B}}^{11}a_{03}+4i {\tilde{B}}^{22}a_{30} \end{aligned}$$

(6.22)

We try to solve the evolution equations analytically in the limit $2 \pi k \omega \gg 1$,along with an additional condition on the natural frequency of the two identical atoms i.e restricting $\omega _o$ to take values such that $\coth (\pi k \omega _o)$=0. The first condition ensures that the factor $\left( 1-e^{-2 \pi k \omega }\right) ^{-1}$ appearing in ${\tilde{A}}^{\alpha \beta }$ and ${\tilde{B}}^{\alpha \beta }$ term of the GSKL matrix $C_{ij}^{\alpha \beta }$ reduces to unity. Substituting the components $H_{ij}^{\alpha \beta }$ and $C_{ij}^{\alpha \beta }$ (calculated in Appendix B and C ) in the evolution equations we obtain the following simplified equations .

$$\begin{aligned} \begin{aligned} {\dot{a}}_{03}(\tau )&= 4B_{1}+A_{1}a_{++}+B_2a_{+-}+B_2a_{-+}-A_1a_{--} \\ {\dot{a}}_{30}(\tau )&= 4B_1+A_1a_{++}+B_2a_{+-}+B_2a_{-+}-A_1a_{--} \\ {\dot{a}}_{++}(\tau )&= 4A_1a_{03}+4A_1a_{30}+2\omega a_{++} \\ {\dot{a}}_{+-}(\tau )&= D_2a_{03}-D_2a_{30}-4B_2a_{30} \\ {\dot{a}}_{-+}(\tau )&= D_2a_{03}-D_2a_{30}-4B_2a_{03}-4B_2a_{30} \\ {\dot{a}}_{--}(\tau )&= -4A_1a_{03}-4A_1a_{30}+2\omega a_{--} \\ {\dot{a}}_{33}(\tau )&= 4B_1 a_{03}+ 4B_2 a_{30} \end{aligned} \end{aligned}$$

(6.23)

where in the above sets of equations the following notations have been used, which are also calculated in Appendix C:

$$\begin{aligned} \begin{aligned} A_1(\omega ) =&\frac{1}{4}(A^{12}+A^{21}) \\ =&-\frac{i \mu ^2}{4L\sqrt{1+(\frac{L}{2k})^2}}~ \cosh \left( (2n+1) ~\mathrm \sinh ^{-1}\left( \frac{L}{2k}\right) \right) \\ B_1(\omega ) =&i {\tilde{B}}^{11}=i \tilde{B^{22}} \\ =&-\frac{\mu ^2}{4\pi k}\left( n+\frac{1}{2}\right) \\ B_2(\omega ) =&i {\tilde{B}}^{12}=i \tilde{B^{21}} \\ =&- \frac{ \mu ^2}{4L\sqrt{1+\left( \frac{L}{2k}\right) ^2}}~ \sinh \left( (2n+1) ~\mathrm \sinh ^{-1}\left( \frac{L}{2k}\right) \right) \\ D_2(\omega ) =&i(B^{12}-B^{21}) = 0 \end{aligned}\nonumber \\ \end{aligned}$$

(6.24)

Solving the Eq. (6.23) we obtain the finite time dependent solution of the Bloch vector components:

$$\begin{aligned} a_{03}(\tau )= & {} C_{1}e^{f_1(\omega ) \tau } \frac{f_1(\omega )}{4(B_1+B_2)} \nonumber \\&+C_2e^{f_2(\omega ) \tau } \frac{f_2(\omega )}{4(B_1+B_2)}\nonumber \\&+C_{3}e^{f_3(\omega ) \tau } \frac{f_3(\omega )}{4(B_1+B_2)} \end{aligned}$$

(6.25)

$$\begin{aligned} a_{30}(\tau )= & {} C_{1}e^{f_1(\omega ) \tau } \frac{f_1(\omega )}{4(B_1+B_2)}\nonumber \\&+C_2e^{f_2(\omega ) \tau } \frac{f_2(\omega )}{4(B_1+B_2)}\nonumber \\&+C_{3}e^{f_3(\omega ) \tau } \frac{f_3(\omega )}{4(B_1+B_2)} \end{aligned}$$

(6.26)

$$\begin{aligned} a_{++}(\tau )= & {} C_{4}e^{2\omega \tau }+C_{1}e^{f_1(\omega ) \tau }\nonumber \\&\left( \frac{-2 f_1(\omega ) A_1}{(-2\omega +f_1)(B_1+B_2)} +\frac{(f_1^2(\omega )+12 B_2^2}{4A_1(B_1+B_2)}\right) \nonumber \\&+e^{f_2(\omega ) \tau }\nonumber \\&\left( \frac{-2 f_2(\omega ) A_1}{(-2\omega +f_2)(B_1+B_2)} +\frac{(f_2^2(\omega )+12 B_2^2}{4A_1(B_1+B_2)}\right) \nonumber \\&+e^{f_3(\omega ) \tau }\nonumber \\&\left( \frac{-2 f_3(\omega ) A_1}{(-2\omega +f_3)(B_1+B_2)} +\frac{(f_3^2(\omega )+12 B_2^2)}{4A_1(B_1+B_2)}\right) \end{aligned}$$

(6.27)

$$\begin{aligned} a_{+-}(\tau )= & {} -C_5-\left( \frac{B_2}{B_1+B_2}\right) \nonumber \\&(C_1 e^{f_1(\omega )\tau }+C_2 e^{f_2(\omega )\tau }+C_3 e^{f_3(\omega )\tau }) \end{aligned}$$

(6.28)

$$\begin{aligned} a_{-+}(\tau )= & {} C_5-\left( \frac{2B_2}{B_1+B_2}\right) \nonumber \\&(C_1 e^{f_1(\omega )\tau }+C_2 e^{f_2(\omega )\tau }+C_3 e^{f_3(\omega )\tau }) \end{aligned}$$

(6.29)

$$\begin{aligned} a_{--}(\tau )= & {} C_4e^{2 \omega \tau }-\frac{2A_1}{B_1+B_2}\nonumber \\&\left( \frac{C_1 f_1(\omega )e^{f_1(\omega )\tau }}{-2\omega +f_1}+ \frac{C_2 f_2(\omega )e^{f_2(\omega )\tau }}{-2\omega +f_2}+\frac{C_3 f_3(\omega )e^{f_3(\omega )\tau }}{-2\omega +f_3}\right) \nonumber \\ \end{aligned}$$

(6.30)

$$\begin{aligned} a_{33}(\tau )= & {} \left( C_6+C_1e^{f_1(\omega )\tau } +C_2e^{f_2(\omega )\tau }+C_3e^{f_3(\omega )\tau }\right) \end{aligned}$$

(6.31)

where $\ C_i~~ \forall i=1,6$ are arbitrary constants.

In the above sets of equations $f_1(\omega )$, $f_2(\omega )$, $f_3(\omega )$ can explicitly be written as:

$$\begin{aligned} \begin{aligned} f_1(\omega )&= -2\omega (\Delta _1(\omega ))^2+(\Delta _1(\omega ))^3-24\omega B_2^2 \\&\quad +\Delta _1(\omega )(-16A_1^2(\omega )+12B_2^2(\omega )) \\ f_2(\omega )&= -2\omega (\Delta _2(\omega ))^2+(\Delta _2(\omega ))^3-24\omega B_2^2\\&\quad +\Delta _2(\omega )(-16A_1^2(\omega )+12B_2^2(\omega )) \\ f_3(\omega )&= -2\omega (\Delta _3(\omega ))^2+(\Delta _1(\omega ))^3-24\omega B_2^2\\&\quad +\Delta _3(\omega )(-16A_1^2(\omega )+12B_2^2(\omega )) \end{aligned} \end{aligned}$$

(6.32)

The functions $\Delta _1(\omega )$,$\Delta _2(\omega )$ and $\Delta _3(\omega )$ appearing in the $f_1(\omega )$,$f_2(\omega )$ and $f_3(\omega )$ are explicitly written in the following equations.

$$\begin{aligned} \begin{aligned} \Delta _1(\omega )&= \frac{2\omega }{3}-2^{1/3} \frac{b_1}{3} {{\mathcal {Z}}}(b_1,b_2) +\frac{1}{3 \times 2^{1/3}} {{\mathcal {Z}}}(b_1,b_2) \\ \Delta _2(\omega )&= \frac{2\omega }{3}+\frac{((1+i \sqrt{3}b_1))}{3\times 2^{2/3} {{\mathcal {Z}}}(b_1,b_2)}\\&\quad -\frac{1}{6\times 2^{1/3}}(1-i \sqrt{3}) {{\mathcal {Z}}}(b_1,b_2) \\ \Delta _3(\omega )&= \frac{2\omega }{3}+\frac{((1-i \sqrt{3}b_1))}{3\times 2^{2/3} {{\mathcal {Z}}}(b_1,b_2)}\\&\quad -\frac{1}{6\times 2^{1/3}}(1+i \sqrt{3}) {{\mathcal {Z}}}(b_1,b_2) \end{aligned} \end{aligned}$$

(6.33)

where we have denoted the following symbols:     

$$ b_1=-48 A_1^2+36B_2^2-4\omega ^2 ~~~~~\ \ \ b_2=288A_1^2\omega +432 B_2^2 \omega +16\omega ^3 $$

and we define again a new function, ${{\mathcal {Z}}} (b_1,b_2) = \left( b_2+\sqrt{4b_1^3+b_2^2}\right) ^{1/3}$. The above equation shows that the function $\Delta _1(\omega )$ is real, whereas $\Delta _2(\omega )$ and $\Delta _3(\omega )$ are complex. The arbitrary constants are determined using the equilibrium behaviour as the boundary conditions which are already mentioned in Eq. 6.10. From the appearance it might seem that the function $e^{f_l(\omega )}$ diverges as $\tau$ tends to infinity but it can be explicitly shown that the function $f_{l}(\omega )(l=1,2,3) <0$, taking into account the leading order term in $\Delta$, which is a decaying function. Hence the function, $e^{f_l(\omega ) \tau }$ tends to a finite value, which we denote by $f_l(\omega )\tau '$.

The physically acceptable solutions to the simplified evolution equations obeying the boundary conditions are given by the following expressions:

$$\begin{aligned} a_{03}(\tau )= & {} -\left[ g_1(\omega )e^{-|f_1(\omega )|(\tau -\tau ')} \frac{|f_1(\omega )|}{4(B_1+B_2)} \right. \nonumber \\&\quad +g_2(\omega )e^{-|f_2(\omega )|(\tau -\tau ')} \frac{|f_2(\omega )|}{4(B_1+B_2)} \nonumber \\&\quad \left. +g_{3}(\omega )e^{-|f_3(\omega )|(\tau -\tau ')} \frac{|f_3(\omega )|}{4(B_1+B_2)}\right] \end{aligned}$$

(6.34)

$$\begin{aligned} a_{30}(\tau )= & {} -\left[ g_{1}(\omega )e^{-|f_1(\omega )| (\tau -\tau ')} \frac{|f_1(\omega )|}{4(B_1+B_2)}\right. \nonumber \\&\quad +g_2(\omega )e^{-|f_2(\omega )| (\tau -\tau ')} \frac{|f_2(\omega )|}{4(B_1+B_2)} \nonumber \\&\quad \left. +g_{3}(\omega )e^{-|f_3(\omega )|(\tau -\tau ')} \frac{|f_3(\omega )|}{4(B_1+B_2)}\right] \end{aligned}$$

(6.35)

$$\begin{aligned} a_{++}(\tau )= & {} g_{1}(\omega )e^{-|f_1(\omega )|(\tau -\tau ')} \nonumber \\&\left( \frac{2|f_1(\omega )| A_1}{-(2\omega +|f_1|)(B_1+B_2)} +\frac{f_1^2(\omega )+12 B_2^2}{4A_1(B_1+B_2)}\right) \nonumber \\&+g_2(\omega ) e^{-|f_2(\omega )| (\tau -\tau ')}\nonumber \\&\left( \frac{2|f_2(\omega )| A_1}{-(2\omega +|f_2|)(B_1+B_2)} +\frac{f_2^2(\omega )+12 B_2^2}{4A_1(B_1+B_2)}\right) \nonumber \\&+g_3(\omega ) e^{-|f_3(\omega )| (\tau -\tau ')}\nonumber \\&\left( \frac{2 |f_3(\omega )| A_1}{-(2\omega +|f_3|)(B_1+B_2)} +\frac{f_3^2(\omega )+12 B_2^2}{4A_1(B_1+B_2)}\right) \nonumber \\ \end{aligned}$$

(6.36)

$$\begin{aligned} a_{+-}(\tau )= & {} -g_5(\omega )-\left( \frac{B_2}{B_1+B_2}\right) \nonumber \\&\left( g_1(\omega ) e^{-|f_1(\omega )|(\tau -\tau ')}+g_2(\omega ) e^{-|f_2(\omega )|(\tau -\tau ')}\right. \nonumber \\&\quad \left. +g_3(\omega ) e^{-|f_3(\omega )|(\tau -\tau ')}|\right) \end{aligned}$$

(6.37)

$$\begin{aligned} a_{-+}(\tau )= & {} g_5(\omega )-\left( \frac{2B_2}{B_1+B_2}\right) \nonumber \\&\left( g_1(\omega ) e^{-|f_1(\omega )|(\tau -\tau ')}+g_2(\omega ) e^{-|f_2(\omega )|(\tau -\tau ')}\right. \nonumber \\&\quad \left. +g_3(\omega ) e^{-|f_3(\omega )|(\tau -\tau ')}\right) \end{aligned}$$

(6.38)

$$\begin{aligned} a_{--}(\tau )= & {} -\frac{2A_1}{B_1+B_2}\nonumber \\&\left[ \frac{g_1(\omega ) |f_1(\omega )|e^{-|f_1(\omega )|(\tau -\tau ')}}{2\omega +|f_1|}\right. \nonumber \\&\quad +\frac{g_2(\omega ) |f_2(\omega )|e^{-|f_2(\omega )|(\tau -\tau ')}}{2\omega +|f_2|} \nonumber \\&\quad \left. +\frac{g_3(\omega ) |f_3(\omega )|e^{-|f_3(\omega )|(\tau -\tau ')}}{2\omega +|f_3|}\right] \end{aligned}$$

(6.39)

$$\begin{aligned} a_{33}(\tau )= & {} g_6(\omega )+g_1(\omega )e^{-|f_1(\omega )|(\tau -\tau ')}\nonumber \\&+g_2(\omega )e^{-|f_2(\omega )|(\tau -\tau ')}+g_3(\omega )e^{-|f_3(\omega )|(\tau -\tau ')}\nonumber \\ \end{aligned}$$

(6.40)

where the explicit functional forms of $g_i(\omega )~\forall i=1,\cdots , 6$ are given below and the constant $C_4$ is zero which is consistent with the given boundary condition.

The arbitrary constants obtained after using the late time behaviour as the initial conditions are explicitly defined by the following expressions:⁵

$$\begin{aligned} g_1(\omega )= & {} \frac{1}{\left( -\frac{\mathcal {Y}_1 f_2(\omega )}{4(B_1+B_2)} +\frac{{\mathcal {Y}}_1 f_3(\omega )}{4(B_1+B_2)}\right) \left( -\frac{{\mathcal {Y}}_2 f_3(\omega )}{{\mathcal {Y}}_4}+\frac{{\mathcal {Y}}_3 f_2(\omega )}{{\mathcal {Y}}_4}\right) +\left( \frac{{\mathcal {Y}}_1 f_2(\omega )}{4(B_1+B_2)^2}+\frac{3B_2 f_3(\omega )}{4(B_1+B_2)^2}\right) \left( -\frac{{\mathcal {Y}}_2 f_3(\omega )}{{\mathcal {Y}}_4}+\frac{{\mathcal {Y}}_3 f_2(\omega )}{{\mathcal {Y}}_4}\right) }\nonumber \\&\times \left[ \frac{-\frac{{\mathcal {Y}}_3}{4A_1(B_1+B_2)}\left( -\frac{{\mathcal {Y}}_1 f_2(\omega )}{4(B_1+B_2)}+\frac{{\mathcal {Y}}_1 f_3(\omega )}{4(B_1+B_2)}\right) \mathrm ~tanh(\pi k \omega )-{\mathcal {Y}}_1(-\frac{{\mathcal {Y}}_2 f_3}{{\mathcal {Y}}_4}+\frac{{\mathcal {Y}}_3 f_2}{{\mathcal {Y}}_4})\mathrm ~tanh(\pi k \omega )}{\left( -\frac{{\mathcal {Y}}_1f_2(\omega )}{4(B_1+B_2)}+\frac{{\mathcal {Y}}_1f_2(\omega )f_3(\omega )}{4(B_1+B_2)}\right) }\right] , \end{aligned}$$

(6.41)

$$\begin{aligned} g_2(\omega )= & {} \frac{4(B_1+B_2)12B_1B_2^2+36B_2^2+B_1f_1^2+3B_2f_3^2 \mathrm ~tanh(\pi k \omega )}{(f_2-f_3)(12 B_1 B_2^2-36B_2^3-B_1f_1^2+B_1f_1f_2+B_1f_1f_3+3B_2f_2f_3)}, \end{aligned}$$

(6.42)

$$\begin{aligned} g_3(\omega )= & {} -\frac{4}{(f_2-f_3)(12 B_1 B_2^2-36B_2^3-B_1f_1^2+B_1f_1f_2+B_1f_1f_3+3B_2f_2f_3)}\nonumber \\&\times [12B_1^2B_2^2 \mathrm{~tanh(\pi k \omega )}+48 B_1B_2^3 \mathrm{~tanh(\pi k \omega )} +36B_2^4\mathrm{~tanh(\pi k \omega )}+B_1^2f_1^2 \mathrm{~tanh(\pi k \omega )}\nonumber \\&+B_1B_2f_1^2\mathrm{~tanh(\pi k \omega )}+3B_1B_2f_2^2\mathrm{~tanh(\pi k \omega )}+3B_2^2f_2\mathrm{~tanh(\pi k \omega )}], \end{aligned}$$

(6.43)

$$\begin{aligned} g_5(\omega )= & {} -\frac{1}{12B_1B_2^2+36B_2^3B_1f_2^2-B_1f_1f_2-B_1f_1f_3-3B_2f_2f_3}\nonumber \\&~~~\times [4B_1B_2f_2~\mathrm{tanh(\pi k\omega )}+3B_2^2f_2~\mathrm{tanh(\pi k\omega )} +B_1B_2f_3~\mathrm{tanh(\pi k \omega )}+3B_2^2~ f_3\mathrm{tanh(\pi k \omega )}], \end{aligned}$$

(6.44)

$$\begin{aligned} g_6(\omega )= & {} \frac{1}{12B_1B_2^2+36B_2^3+B_1f_1^2-B_1f_1f_2-B_1f_1f_3-3B_2f_2f_3}\nonumber \\&\times [4B_1^2f_2~\mathrm{tanh(\pi k\omega )}\mathrm{tanh(\pi k\omega )}+16B_1B_2(f_2~\mathrm{tanh(\pi k\omega )}+12B_2^2f_3~\mathrm{tanh(\pi k\omega )}\nonumber \\&+4B_1^2f_5~\mathrm{tanh(\pi k\omega )}+16B_1B_2~\mathrm{tanh(\pi k\omega )}+12B_2^2f_3~\mathrm{tanh(\pi k\omega )}\nonumber \\&-12B_1B_2^2~\mathrm{tanh^2(\pi k\omega )}-36B_2^3~\mathrm{tanh^2(\pi k\omega )}-B_1f_1^2~\mathrm{tanh^2(\pi k\omega )}\nonumber \\&+B_1f_1f_2~\mathrm{tanh^2(\pi k\omega )}+B_1f_1f_3~\mathrm{tanh^2(\pi k\omega )}\nonumber \\&+3B_2f_2f_3~\mathrm{tanh^2(\pi k\omega )}]. \end{aligned}$$

(6.45)

where in writing the function $g_1(\omega )$ we have introduced the following symbols:

$$\begin{aligned} {\mathcal {Y}}_1= & {} \left( 1-\frac{B_2}{B_1+B_2}\right) ,~~ {\mathcal {Y}}_2=-12B_2^2-f_2^2, \nonumber \\ {\mathcal {Y}}_3= & {} -12B_2^2-f_3^2,~~ {\mathcal {Y}}_4=16A_1(B_1+B_2)^2. \end{aligned}$$

(6.46)

Using these solution the variation of the Bloch vector components with respect to various parameters are plotted in Figs. 3, 4, 5 and 6.

Fig. 3

Dependence of the Bloch vector components on Time is shown here

Fig. 4

Dependence of the Bloch vector components on frequency is shown here

Fig. 5

Dependence of the various Bloch vector components on the Euclidean distance is shown here

Fig. 6

Dependence of the Bloch vector components on k is shown here

7 Entanglement measures

Using these solution our next objective is to compute various entanglement measures from the present OQS set up.

7.1 Von Neumann entanglement entropy

In the context of quantum statistical mechanics and quantum information theory, Von Neumann entropy plays the role of extended version of classical Gibbs entropy. It actually measures the amount of quantum entanglement for a subsystem or reduced system of a bipartite quantum system.

In the present context, for the two atomic subsystem, which we have obtained after partially tracing over the bath degrees of freedom the Von Neumann entanglement entropy is defined in terms of the reduced density matrix ($\rho _\mathbf{System }$) as:

$$\begin{aligned} S(\rho _\mathbf{System })=-\mathrm{Tr}\left[ \rho _\mathbf{System }~\ln \left( \rho _\mathbf{System }\right) \right] \end{aligned}$$

(7.1)

• 1. It is expected from our present OQS set up that the measure of Von Neumann entanglement entropy is non zero because our subsystem, which we have obtained by partially tracing over bath degrees of freedom, is described by mixed states. Conversely, if a subsystem is to be characterised by pure quantum mechanical states , the entanglement measure is zero. Since in our set up we are dealing with pure quantum mechanical states initially due to the initial uncorrelation between the system and the environment, we get zero entanglement measure from our computation. But during time evolution of the subsystem it becomes more correlated and consequently pure state transforms to a mixed state, for which we get a saturation but non-zero value of the Von Neumann entanglement entropy.

• 2. The Von Neumann entanglement measure for the reduced subsystem can be expressed as a sum of the contributions from the two independent and identical atoms.

  $$\begin{aligned} S(\rho _\mathbf{System }) = S(\rho _1 \otimes \rho _2) = S(\rho _1)+S(\rho _2) \end{aligned}$$

  (7.2)

• 3. Additionally, it is important to note that, the Von Neumann entanglement measure is not the only reasonable measure of quantum entanglement. In the later sections we will consider other quantum entanglement measures, which are commonly used in the context of quantum information theory.

In the context of two entangled atoms relevant to our OQS it can be seen from the solution of the Bloch vector components that $\mathrm {a_{03}}(\tau )$ and $\mathrm {a_{30}}(\tau )$ has the same form i.e. $a_{03}(\tau )=a_{30}(\tau )$. So in that case the Von Neumann entropy can be written as:

$$\begin{aligned} S= & {} \frac{1}{4} [ 8~\ln ~2-(\alpha -{\widetilde{\beta }}) \ln (\alpha -{\widetilde{\beta }}) - (\alpha +{\widetilde{\beta }}) \ln (\alpha +{\widetilde{\beta }}) \nonumber \\&\quad - (\eta -{\widetilde{\gamma }}) \ln (\eta -{\widetilde{\gamma }}) - (\eta +{\widetilde{\gamma }}) \ln (\eta +{\widetilde{\gamma }}) ] \end{aligned}$$

(7.3)

where we have introduced two new functions, ${\widetilde{\beta }}$ and ${\widetilde{\gamma }}$, which are defined as:

$$\begin{aligned} {\widetilde{\beta }}= & {} \beta _\mathbf{Two atom }=\sqrt{a_{-+}(\tau )a_{+-}(\tau )}, \nonumber \\ {\widetilde{\gamma }}= & {} \gamma _\mathbf{Two atom }=\sqrt{4a^2_{03}(\tau )+ a_{--}(\tau )a_{++}(\tau )} \end{aligned}$$

(7.4)

Here all of these time dependent coefficients $a_{03}(\tau )$, $a_{-+}(\tau )$, $a_{+-}(\tau )$, $a_{--}(\tau )$, $a_{++}(\tau )$ and $a_{33}(\tau )$ have explicitly been computed for the two atomic OQS in the previous section.

Fig. 7

Normalized Von Neumann Entanglement entropy variation with various parameters is shown here

To physically analyse this result we have plotted the behaviour of the Von Neumann entropy from the present two atomic OQS set up with respect to different useful parameters present in the theory. However, as we will see later that almost all the netanglement measures give plots of similar nature. Hence, the discussions related to all the plots are given afterwards in a separate Sect. 9.

7.2 R$e^{'}$nyi entropy

R$e^{'}$nyi entropy is a generalisation of various information theoretic measure which quantify randomness of a quantum mechanical system. The R$e^{'}$nyi entropy can be expressed in terms of Hartley function as:

$$\begin{aligned}&S_q(\rho _\mathbf{System })=\frac{1}{1-q}\ln \mathrm{Tr} [(\rho _\mathbf{System })^q]\equiv {\mathcal {H}}_{q}(\rho _\mathbf{System }), \nonumber \\&\quad q\ge 0,q\ne 1 \end{aligned}$$

(7.5)

where q is known as the R$e'$nyi Index.

In the case of our two atomic OQS set up, we found from the structure of the solution that $a_{03}$=$a_{30}$. Hence the expression for the R$e'$nyi entropy reduces to the following simplified form:

$$\begin{aligned}&S_q(\rho _\mathbf{System }) \nonumber \\&\quad = \frac{1}{1-q} \ln \nonumber \\&\qquad \left[ 4^{-q} \left\{ (\alpha -{\widetilde{\beta }})^q + (\alpha +{\widetilde{\beta }})^q + (\eta -{\widetilde{\gamma }})^q + (\eta +{\widetilde{\gamma }})^q \right\} \right] \nonumber \\ \end{aligned}$$

(7.6)

where the symbols used has already been defined in the previous section. It can be very easily verified that in this context relevant to our OQS the R$e'$nyi entropy reduces to the Von Neumann entropy in the limit q $\rightarrow 1$.

Fig. 8

Renyi entropy (q $\rightarrow$1) variation with various parameters are shown here

Fig. 9

Collision entropy variation with various parameters are shown here

Fig. 10

Min entropy variation with various parameters are shown here

7.3 Logarithmic negativity

Let us start with Peres–Horodeski criterion which is the necessary condition for the joint density matrix of two quantum mechanical system A and B to be separable. This is sometimes called Positive Partial Transpose (PPT) criterion. This is mainly used to decide the separability of mixed states, where Schmidt decomposition does not apply. Using this criteria one can defineLogarithmic Negativity as:

$$\begin{aligned} E_{N}(\rho _\mathbf{System } )=\ln ||\rho ^{T_A}|| \end{aligned}$$

(7.7)

where, $||\rho ^T_A||$ is the trace norm and is defined as:

$$\begin{aligned} ||\rho ^T_A||= & {} \mathrm{Tr} \left( \sqrt{ (\rho ^T_A)^\dagger (\rho ^T_A)}\right) =\sum _i |\lambda _i| \nonumber \\= & {} \sum _{\lambda _i>0} \lambda _i+\sum _{\lambda _i<0}|\lambda _i|=2\sum _{\lambda _i<0}|\lambda _i|+1=2N+1.\nonumber \\ \end{aligned}$$

(7.8)

Then the logarithmic negativity can be recast as:

$$\begin{aligned} E_{N}(\rho _\mathbf{System } )=\ln (2N+1). \end{aligned}$$

(7.9)

If we fix $N=0$ then it represents no entanglement in the quantum system. For this non entangled case one can write the system density matrix as:

$$\begin{aligned} \rho _\mathbf{System } =\sum _i \lambda _i \rho _i^A \otimes \rho _i^B=\rho ^{T_A}, \end{aligned}$$

(7.10)

where we define the sub system density matrix as:

$$\begin{aligned} \rho _i^q=|i\rangle _{q}{}_{q}\langle i| ~~\forall ~~q=A,B~~~\mathrm{with}~~ \lambda _i \ge 0 \end{aligned}$$

(7.11)

For $N\ne 0$ which represents the entangled case the system density matrix can be expressed as:

$$\begin{aligned} \rho _\mathbf{System }= & {} \sum _{i,j,k,l} {\mathcal {C}}_{ijkl}\left( |i\rangle _{AA}\langle j|\right) \otimes \left( |k\rangle _{AA}\langle l|\right) \nonumber \\= & {} \sum _{i,j,k,l} {\mathcal {C}}_{ijkl}(|j\rangle _{AA}\langle i|) \otimes (|k\rangle _{AA}\langle l|) \end{aligned}$$

(7.12)

So from this computation if we get negative eigen values for a quantum mechanical system then we can say that the corresponding quantum states are entangled.

For the case of two entangled atoms, using the solutions of the bloch vector components, the logarithmic negativity is given by

$$\begin{aligned}&E_N(\rho )=\ln \left[ \frac{17}{100}\left\{ (2+8a^2_{03})+a_{33}(4+2a_{33})+a^2_{-+}+a^2_{+-}\right. \right. \nonumber \\&\quad \left. \left. \times -\sqrt{(4+10a_{33}+(a_{-+}-a_{+-})^2)(16a^2_{03})+(a_{-+}+a_{+-})^2)}\right\} ^{1/2} \right. \nonumber \\&\quad \left. +\frac{17}{100}\left\{ (2+8^2a_{03})+a_{33}(4+2a_{33})+a^2_{-+}+a^2_{+-}\right. \right. \nonumber \\&\quad \quad \left. \left. \times +\sqrt{(4+10a_{33}+(a_{-+}-a_{+-})^2)(16a^2_{03}+(a_{-+}+a_{+-})^2)}\right\} ^{1/2} \right. \nonumber \\&\quad \left. +\frac{17}{100}\left\{ (2+a_{33}(-4+2a_{33})+a^2_{--}+a^2_{++}\right. \right. \nonumber \\&\quad \left. \left. \times -\sqrt{(4+6a_{33}+(a_{--}-a_{++})^2)((a_{--}+a_{++})^2))}\right\} ^{1/2} \right. \nonumber \\&\quad \left. +\frac{17}{100}\left\{ (2+a_{33}(-4+2a_{33})+a^2_{--}+a^2_{++} \right. \right. \nonumber \\&\quad \left. \left. +\sqrt{(4+6a_{33}+(a_{--}-a_{++})^2)(4(\beta ^2-a_{-+}a_{+-}+(a_{--}+a_{++})^2))}\right\} ^{1/2}\right] \nonumber \\ \end{aligned}$$

(7.13)

7.4 Entanglement of formation and concurrence

Both of the measures studied in earlier sections are used to quantify the resources needed to create a given entangled state. Each of them are used as an entanglement measure for bipartite quantum state in quantum information theory. The entanglement of formation for pure and mixed states takes the following form:

$$\begin{aligned} \begin{array}{lll} \displaystyle E_{f}(\rho ) =\left\{ \begin{array}{ll} \displaystyle -\mathrm{Tr}(\rho _{A}\ln \rho _{A})=- \mathrm{Tr}(\rho _{B}\ln \rho _{B})~~~~ &{} {\mathbf{for pure state }} \\ \displaystyle \mathbf{inf}\left( \sum _j p_jE_f(\Phi _j)\right) ~~~~ &{} {\mathbf{for mixed state }} \\ \end{array} \right. \end{array} \end{aligned}$$

(7.14)

For mixed states the infimum is taken over all possible decompositions of the density matrix $\rho$ into pure states. The quantities used in the above equation are defined below:

$$\begin{aligned} \rho _A= & {} Tr_B \rho ,~~~ \rho _B= Tr_A \rho , \nonumber \\ \rho _\mathbf{System }= & {} |\Phi \rangle \langle \Phi |~~~~~~~ {\mathbf{for pure state }} \end{aligned}$$

(7.15)

$$\begin{aligned} E_f(\Phi _j)= & {} -\mathrm{Tr}(\Phi _jln \Phi _j), \nonumber \\ \rho _\mathbf{System }= & {} \sum _j p_j |\Phi _j\rangle \langle \Phi _j|~~~~~~ {\mathbf{for mixed state }} \end{aligned}$$

(7.16)

Relation between entanglement of formation and Concurrence can be written as:

$$\begin{aligned} E_f(\rho _\mathbf{System })= & {} {\mathcal {E}}(C(\rho _\mathbf{System })) \nonumber \\= & {} h\left( \frac{1+\sqrt{1-C^2(\rho _\mathbf{System })}}{2}\right) \end{aligned}$$

(7.17)

where h(x) is known as the Binary Entropy function which is defined as:

$$\begin{aligned} h(x)=-x\ln x-(1-x)\ln (1-x) \end{aligned}$$

(7.18)

Fig. 11

Renyi entropy variation with Renyi Index is shown here

Fig. 12

Log Negativity variation with various parameters are shown here

Fig. 13

Variation of the binary entropy function h(x) with x is shown here

Here we study about concurrence in the context of two entangled atoms from the perspective of OQS. For X type states it is defined analytically by the following expression:

$$\begin{aligned} C(\rho _\mathbf{System })=\mathrm{max}[0,\lambda _1-\lambda _2-\lambda _3-\lambda _4] \end{aligned}$$

(7.19)

where the $\lambda _i$’s are the square roots of the eigenvalues of the matrix:

$$\begin{aligned} \begin{array}{lll} \displaystyle {\mathcal {R}}\equiv \left\{ \begin{array}{ll} \displaystyle \sqrt{\sqrt{\rho _\mathbf{System }}~{\tilde{\rho }}_\mathbf{System }~\sqrt{\rho _\mathbf{System }}}~~~~ &{} {\mathbf{for Hermitian }} \\ \displaystyle \sqrt{\rho _\mathbf{System } ~{\tilde{\rho }}_\mathbf{System }}~~~~ &{} {\mathbf{for Non Hermitian }} \\ \end{array} \right. \end{array} \end{aligned}$$

(7.20)

where ${\tilde{\rho }}_\mathbf{System }$ is the spin flip (Werner type) quantum states [48] given by⁶:

$$\begin{aligned} {\tilde{\rho }}_\mathbf{System }= & {} (\sigma _2 \otimes \sigma _2) \rho ^*(\sigma _2 \otimes \sigma _2) \nonumber \\= & {} [((\sigma _--\sigma _+)\otimes (\sigma _--\sigma _+)\rho ^*(\sigma _-- \sigma _+)\otimes (\sigma _--\sigma _+))]\nonumber \\ \end{aligned}$$

(7.23)

and $\rho ^*_\mathbf{System }$ in the above expression indicates complex conjugate of $\rho _\mathbf{System }$.

The eigenvalues $\lambda$’s follow the following sequence:

$$\begin{aligned} \lambda _1>\lambda _2>\lambda _3>\lambda _4, \mathrm{where}~~~ \lambda _1-\lambda _2-\lambda _3-\lambda _4>0 \end{aligned}$$

(7.24)

Fig. 14

Relation between the entanglement of formation and concurrence is shown here

For the case of two entangled atoms relevant to our system the solutions of the Bloch vector components ${ a_{30}}$ and ${ a_{03}}$ obtained are identical. Thus the eigenvalues $\lambda _1$, $\lambda _2$, $\lambda _3$ and $\lambda _4$ are given by the following expressions:

$$\begin{aligned} \lambda _1= & {} \frac{1}{4}\sqrt{{{{\mathcal {S}}}}-2{{\mathcal {F}}}} ,~~~ \lambda _2=\frac{1}{4}\sqrt{{{\mathcal {S}}}+2{{\mathcal {F}}}}, \nonumber \\ \lambda _3= & {} \frac{1}{4}\sqrt{{{\mathcal {X}}}-2{{\mathcal {U}}}}, ~~~ \lambda _4=\frac{1}{4}\sqrt{{{\mathcal {X}}}+2{{\mathcal {U}}}}.~~~ \end{aligned}$$

(7.21)

where we define ${{\mathcal {S}}}$, ${{\mathcal {X}}}$, ${{\mathcal {F}}}$ and ${{\mathcal {U}}}$ as:

$$\begin{aligned} \begin{aligned} {{\mathcal {S}}}&= 1-2a_{33}+a^2_{33}+a_{-+}a_{+-} \\ {{\mathcal {X}}}&=1-4a^2_{03}+2a_{33}+a^2_{33}+a_{--}a_{++} \\ {{\mathcal {F}}}&= \sqrt{a_{-+}a_{+-}-2a_{33}a_{-+}a_{+-}+a^2_{33}a_{-+}a_{+-}} \\ {{\mathcal {U}}}&= \sqrt{a_{--}a_{++}-4a^2_{03}a_{--}a_{++} +2a_{33}a_{--}a_{++}+a^2{33}a_{--}a_{++}} \end{aligned} \end{aligned}$$

(7.25)

Therefore, using Eq. 7.19 the concurrence for two entangled atoms relevant to our system can be calculated as:

$$\begin{aligned}&C(\rho _\mathbf{System }) \nonumber \\&\quad =\mathrm{max}\left[ 0,\frac{1}{4}\left( \sqrt{{{\mathcal {S}}}-2{{\mathcal {F}}}}-\sqrt{{\mathcal {S}} +2{{\mathcal {F}}}}-\sqrt{{{\mathcal {X}}}-2{{\mathcal {U}}}}-\sqrt{{{\mathcal {X}}}+ 2{\mathcal {U}}}\right) \right] \nonumber \\ \end{aligned}$$

(7.26)

Fig. 15

Variation of Concurrence with various parameters is shown here

The explicit expression for the entanglement of formation for our system the is given by the following expression:

$$\begin{aligned} E_{f}(\rho )= & {} \frac{1}{2}\left( -1-\sqrt{1-(C(\rho ))^2}\right) \nonumber \\&\log \left[ \frac{1}{2}\left( 1+\sqrt{1-C(\rho )}\right) \right] \nonumber \\&-\left( 1+\frac{1}{2}\left( -1-\sqrt{1-C(\rho )}\right) \right) \nonumber \\&\log \left[ 1+\frac{1}{2}\left( -1-\sqrt{1-C(\rho )}\right) \right] \end{aligned}$$

(7.27)

The following section shows various plots of entanglement of formation calculated from concurrence of our model.

Fig. 16

Variation of entanglement of formation with various parameters is shown here

7.5 Quantum discord

It is a measure of non classical correlation between two subsystems of a quantum system. It includes correlations that are due to quantum physical effects, but do not necessarily involve the concept of quantum entanglement. Sometimes it is also identified as measure of quantumness of correlation functions. If the two quantum states are separable then it does not imply that the quantum correlations exist between them. It is defined as

$$\begin{aligned} {\mathcal {D}}_A(\rho _\mathbf{System })={\mathcal {I}}(\rho _\mathbf{System })-\mathbf{max}_{\Pi _j^A}{\mathcal {T}}_{\Pi _j^A(\rho _\mathbf{System })} \end{aligned}$$

(7.28)

where mutual information ${\mathcal {I}}(\rho _\mathbf{System })$ is defined as:

$$\begin{aligned} {\mathcal {I}}(\rho _\mathbf{System })=S(\rho _A)+S(\rho _B)-S(\rho _\mathbf{System }) \end{aligned}$$

(7.29)

where $S(\rho _A)$, $S(\rho _B)$ and $S(\rho _\mathbf{System })$ represent the Von Neumann entropy of system A (Atom 1), $B (\mathbf{Atom~2})$ and the combined system respectively. On the other hand, the part of the correlation that can be attributed to the classical correlation, which is represented by, $_{\Pi _j^A}{\mathcal {T}}_{\Pi _j^A}(\rho _\mathbf{System })$, is defined as:

$$\begin{aligned} {}_{\Pi _j^A}{\mathcal {T}}_{\Pi _j^A}(\rho _\mathbf{System })= & {} S(\rho _B)-\Pi _j^A S(\rho _B|\Pi _j^A) \end{aligned}$$

(7.30)

Now we know that for any two qubit state the density matrix is given by the following expression:

$$\begin{aligned} \rho =\frac{1}{4}\left( I^a \otimes I^b+\sum _{i=1}^3(a_i\sigma _i \otimes I^b+I^a \otimes b_i\sigma _i)+\sum _{i,j=1}^3 T_{ij}\sigma _i \otimes \sigma _j \right) .\nonumber \\ \end{aligned}$$

(7.31)

Then the geometric measure of quantum discord is evaluated as:

$$\begin{aligned} {{\mathcal {D}}}(\rho )=\frac{1}{4}(||a||^2+||T||^2-\lambda _{max}) \end{aligned}$$

(7.32)

where a is the column vector, which is defined as, $a=(a_1+a_2+a_3)^t$. Here the superscript t stands for the transpose of the vectors or matrices, the trace norm square is defined as:$||a||^2=\sum _i a_i^2$ and $T=(t_{ij})$ is a matrix which one can compute for a specific quantum system and $\lambda _{max}$ is the largest eigenvalue of the matrix $(aa^t+TT^t)$. The matrix T and the column vector relevant to our system are given by:

$$\begin{aligned} T=\begin{pmatrix} a_{++} &{} &{} a_{+-} &{} &{} 0\\ a_{-+} &{} &{} a_{--} &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} a_{33} \\ \end{pmatrix}~~~~~~~~~~ a=\begin{pmatrix} 0\\ 0 \\ a_{30} \\ \end{pmatrix} \end{aligned}$$

(7.33)

Using the above definition, the square of the norm of the matrix T (The norm of any matrix M is given by: ||M||=$\sqrt{\mathrm{Tr}(M^\dagger M)}$) and the column vector a relevant to the system studied is given by:

$$\begin{aligned} ||T||^2= & {} a^2_{33}+a^2_{--}+a^2_{-+}+a^2_{+-}+a^2_{++} \nonumber \\ ||a||^2= & {} a^2_{30}. \end{aligned}$$

(7.34)

The matrix $(\mathrm {aa^t+TT^t})$ then becomes

$$\begin{aligned} (\mathrm {aa^t+TT^t}) =\begin{pmatrix} a^2_{+-}+a^2_{++} &{} &{}a_{--}a_{+-}+a_{-+}+a_{++} &{} &{} 0\\ a_{--}a_{+-}+a_{-+}+a_{++} &{} &{}a^2_{--}+a^2_{-+} &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} a^2_{33}+a^2_{30} \\ \end{pmatrix}\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \end{aligned}$$

(7.35)

with eigenvalues

$$\begin{aligned} \lambda _1=a^2_{33}+a^2_{03}~~~~ \lambda _2=\frac{1}{2}({{\mathcal {P}}}-{\mathcal {Q}}) ~~~~ \lambda _3=\frac{1}{2}({{\mathcal {P}}}+{{\mathcal {Q}}}) \end{aligned}$$

(7.36)

where ${{\mathcal {P}}}$ and ${{\mathcal {Q}}}$ are

$$\begin{aligned} {{\mathcal {P}}}= & {} a^2_{--}+a^2_{-+}+a^2_{+-}+a^2_{++} \end{aligned}$$

(7.37)

$$\begin{aligned} {{\mathcal {Q}}}= & {} \sqrt{(-a^2_{--}-a^2_{-+}-a^2_{+-}-a^2_{++})^2-4(a^2_{-+}a^2_{+-}-2a_{--}a_{-+}a_{+-}a_{++}+a^2_{--}a^2_{++})} \end{aligned}$$

(7.38)

It has already been mentioned in the earlier sections that the solution of the Bloch vector components $a_{30}$ and $a_{03}$ are identical.

Thus the quantum discord for the two entangled atoms relevant to our system calculated using Eq. 7.32 is therefore given by the following expression:

$$\begin{aligned} {\mathcal {D}}(\rho )= & {} \frac{1}{4}\left[ a^2_{03}+a^2_{33}+{\mathcal {Q}}\right] \nonumber \\&-\mathrm{max}\left[ a_{03}^2+a^2_{33},\frac{1}{2}({\mathcal {Q}}- \sqrt{{\mathcal {Q}}^2-4{\mathcal {W}}}), \right. \nonumber \\&\qquad \left. \frac{1}{2}({\mathcal {Q}} +\sqrt{{\mathcal {Q}}^2-4{\mathcal {W}}})\right] \end{aligned}$$

(7.39)

where the symbols ${\mathcal {Q}}$ and ${\mathcal {W}}$ are given as follows:

$$\begin{aligned} {\mathcal {Q}}= & {} a^2_{--}+a^2_{-+}+a^2_{+-}+a^2_{++} \nonumber \\ {\mathcal {W}}= & {} a^2_{-+}a^2_{+-}-2a_{--}a_{-+}a_{+-}a_{++}+a^2_{--}a^2_{++} \end{aligned}$$

(7.40)

For a given set of parameters the maximum eigenvalue is calculated and the following set of plots is obtained by varying various parameters appearing in our model.

Fig. 17

Variation of the geometric measure of Quantum discord with various parameters is shown here

In Fig. 17a, we have explicitly shown the behaviour of quantum discord of our two atomic OQS set up in static patch of de Sitter space with respect to rescaled time T. We have normalized the values with the result obtained from $T=10000$ to properly interpret the obtained result from our model. We consider both small and large time scale limiting situations in this context. It is clearly observed that initially the value of quantum discord is almost zero implying that our two atomic OQS set up do not show any signature of quantum correlations. As time goes on, the sub-system gets more and more quantum mechanically correlated and at a very late time scale, it almost saturates to unity, which implies the maximum measure one can obtain from our model to get quantum correlation. Now as we have passed the test for Von Neumann entropy for our model i.e. that this measure is non-zero then one can surely say that non-zero value of quantum discord and Von Neumann entropy together imply the existence of quantum entanglement in the context of our present model of discussion.

In Fig. 17b, we have shown the behaviour of quantum discord with respect to frequency $\omega$ by keeping all other parameters of the model fixed. We have normalized the values with the result obtained from $\omega =5\times 10^{-2}$. For better understanding the underlying physics we consider both small and large frequency scale limiting situations. It is clearly observed that initially for a finite frequency scale, the value of quantum discord is almost constant showing maximum measure of quantum correlations obtained from quantum discord. From this figure it is easily observed that, as frequency is further increased the obtained measure of quantum discord decreases for our system implying reduction in quantum correlations and when it gets close to $\omega _0$ the value of quantum discord is almost approximately zero indicating no quantum correlation in our system.

In Fig. 17c, we have shown the behaviour of quantum discord with respect to the Euclidean distance L. Again like previous plots, we have normalized the values with the result obtained from $L=10^4$ to explain the underlying physics from this system. We consider both small and large length scale limiting situations. It is observed that initially there is some fluctuation with increase in quantum discord for small length scale. But as the Euclidean distance between the two atoms increase, the subsystem becomes quantum correlated to a maximum saturated value. With further increase in Euclidean distance between the two atoms, there is no change in value implying that the quantum correlation between the atoms in our OQS set up reaches its maximum value.

Last but not the least, in Fig. 17d, we have shown the behaviour of quantum discord with respect to the parameter k, which is basically proportional to the curvature scalar or the Ricci scalar of the static patch of de Sitter space. Like previous plots here also we have normalized the obtained value of quantum discord with the result obtained from $k=10^4$. Similarly like previously mentioned all the plots we consider both small and large values of the k scale limiting situations to interpret the obtained result from our OQS set up. It is clearly observed that initially there are fluctuations in quantum spectrum discord with increase in the parameter k. With further increase in the value of k, the value remains constant implying no change in quantum correlations which reaches its accessible stable very small value. This further implies that the effect of quantum correlation is extremely small for the large value of the parameter k.

8 Non locality from bell CHSH inequality in de Sitter space

The authors have studied Bell violation in Quantum Mechanics and Primordial cosmology in their earlier papers [49,50,51,52,53]. Here, in analogy to those, we generalise to the Bell violation in dS spacetime.

To establish the concept of non-locality in de Sitter space let us start with a quantum mechanical Bell CHSH operator, which can be defined in the following form:

$$\begin{aligned} {{{\mathcal {B}}}}_\mathbf{CHSH } = \left( \mathbf{a}\cdot \sigma \right) \otimes \left\{ (\mathbf{b} + \mathbf{b}').\sigma \right\} + \left( \mathbf{a'}\cdot \sigma \right) \otimes \left\{ (\mathbf{b}-\mathbf{b}') \cdot \sigma \right\} \nonumber \\ \end{aligned}$$

(8.1)

where $\mathbf{a},\mathbf{b},\mathbf{a'}$ and $\mathbf{b'}$ are real unit vectors which play significant role to establish non-locality in the present context.

Now the Bell CHSH inequality states that,⁷ $|\langle B_\mathbf{CHSH }\rangle | \le 2$. To establish the non-locality we necessarily have to violate CHSH inequality in de Sitter spacetime, which is of course not a trivial task to do. The main problem in de Sitter space to generate the effect of long range quantum correlation at late time scale from a non-local Bell’s inequality violating set up. However, in Refs. [49,50,51,52,53] we and other authors have explicitly shown that in case of axion one can construct such a Bell’s inequality violating set up, where the axion effective potential is generated from string theory. In the present set up instead of choosing axion as a Bell’s inequality violating candidate to establish non-locality in quantum mechanics we establish this in a more general and model independent way. We use the density matrix formalism from quantum statistical mechanics where the general density matrix for a quantum system can be parametrized as:

$$\begin{aligned} \rho =\frac{1}{4}\left[ I\otimes I+\mathbf{a \cdot \sigma }\otimes I+I\otimes \mathbf{b.\sigma }+\sum _{j,k=1}^3 c_{jk}\sigma _j \otimes \sigma _k\right] \nonumber \\ \end{aligned}$$

(8.2)

where $\mathbf{a}$, $\mathbf{b}$ and $c_{jk}$ all are in general time dependent quantities which can be explicitly obtained by solving the GSKL master equation in presence of the effective Hamiltonian and the quantum dissipator Lindbladian operator . Here, the vectors $\mathbf{a}$ and $\mathbf{b}$ are given by the following general representation in terms of the elements of the density matrix:

$$\begin{aligned} \mathbf{a}= & {} (0,0,\rho _{11}+\rho _{22}-\rho _{33}-\rho _{44}) \nonumber \\ b= & {} (0,0,\rho _{11}+\rho _{33}-\rho _{22}-\rho _{44}) \end{aligned}$$

(8.3)

and the $c_{jk}$ matrix takes the form

$$\begin{aligned} c_{jk}=\frac{1}{4}\begin{pmatrix} 2\rho _{23}-2(\rho _{14}) &{} &{} 2(\rho _{14})_I &{} &{} 0\\ 2(\rho _{14})_I &{} &{} 2\rho _{23}-2(\rho _{14}) &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} \rho _{11}+\rho _{44}-\rho _{22}-\rho _{33} \\ \end{pmatrix}\nonumber \\ \end{aligned}$$

(8.4)

In the context of OQS described by the two entangled atoms the vectors a and b are given by the simplified form, $\mathbf{a}=(0,0,a_{30})$, $\mathbf{b}=(0,0,a_{03}),$ and the coefficient matrix $c_{jk}$ can be written as:

$$\begin{aligned} c_{jk}=\frac{1}{4}\begin{pmatrix} 2a_{+-}-a_{++}-a_{--} &{} &{} i(a_{--}-a_{++}) &{} &{} 0\\ i(a_{--}-a_{++}) &{} &{} 2a_{+-}+a_{++}+a_{--} &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} a_{33} \\ \end{pmatrix} \end{aligned}$$

(8.5)

Now, Bell CHSH inequality is violated in OQS if and only if when the sum of the two largest eigenvalues of the following matrix satisfy the following constraint condition, $c~ c^\dagger >1$, where the eigenvalues are:

$$\begin{aligned} \lambda _1= & {} 4\left( \rho _{11}+\rho _{44}-\frac{1}{2}\right) ^2=a^2_{33} \nonumber \\ \lambda _{2,3}= & {} 4(|\rho _{14}|\pm \rho _{23})^2=\frac{1}{4}(|a_{--}|\pm a_{+-})^2 \end{aligned}$$

(8.6)

For the initial separable state $\rho _0=|00\rangle \langle 00|$, we cannot expect these eigenvalues to exceed unity after evolution because only non-zero component of the initial state is $\rho _{44}=1$. However, the Bell CHSH inequality provides only a necessary condition for the LHV (local hidden variable) description and does not guarantee existence of an LHV. For this reason we need to pass each detector through a local filter, which transforms the matrix c and $\rho$ in the following form.

$$\begin{aligned} c'= & {} (f_{A}\otimes f_{B})~c~(f_{A}\otimes f_{B}) \nonumber \\ \rho '= & {} (f_{A}\otimes f_{B})~\rho ~(f_{A}\otimes f_{B}) \end{aligned}$$

(8.7)

where the two local filters $f_A$ and $f_B$ are described by the following square matrix:

$$\begin{aligned} f_{A}=f_{B}=\begin{pmatrix} 1 &{} &{} &{} &{} &{} &{} 0\\ 0 &{} &{} &{} &{} &{} &{} \eta \\ \end{pmatrix} \end{aligned}$$

(8.8)

The matrices $\rho '$ and $c'$ can be represented as

$$\begin{aligned} \rho '= & {} \frac{1}{\rho _{11}+\eta ^2(\rho _{22}+\rho _{33})+\eta ^4 \rho _{44}} \nonumber \\&\begin{pmatrix} \rho _{11} &{} &{} 0 &{} &{} 0 &{} &{} \eta ^2 \rho _{14}\\ 0 &{} &{} \eta ^2 \rho _{22} &{} &{} \eta ^2 \rho _{22} &{} &{} 0 \\ 0 &{} &{} \eta ^2 \rho _{23} &{} &{} \eta ^2 \rho _{33} &{} &{} 0 \\ \eta ^2 \rho ^*_{14} &{} &{} 0 &{} &{} 0 &{} &{} \eta ^4 \rho _{44} \end{pmatrix}, \end{aligned}$$

(8.9)

$$\begin{aligned} c'= & {} \frac{2\eta ^2}{\rho _{11}+\eta ^2(\rho _{22}+\rho _{33})+\eta ^4 \rho _{44}} \nonumber \\&\begin{pmatrix} \rho _{23}-(\rho _{14})_R &{} &{} (\rho _{14})_R &{} &{} 0 \\ (\rho _{14})_R &{} &{} \rho _{23}+(\rho _{14})_R &{} &{} 0 \\ 0 &{} &{} 0 &{} &{} \frac{\rho _{11}-\eta ^2(\rho _{22}+\rho _{33})+\eta ^4 \rho _{44}}{2\eta ^2} \end{pmatrix}.\nonumber \\ \end{aligned}$$

(8.10)

Fig. 18

Verification of the Violation of the Bell-CHSH inequality with various parameters is shown here

The eigenvalues of the new matrix $c'(c')^{\dagger }$ are appended below:

$$\begin{aligned} \begin{aligned} \lambda _1'&=\frac{\rho _{11}-\eta ^2(\rho _{22}+\rho _{33})+\eta ^4 \rho _{44}}{\rho _{11}+\eta ^2(\rho _{22} +\rho _{33})+\eta ^4 \rho _{44}} \\&=\frac{(1-2\eta ^2)+(1-\eta ^4)(a_{03}+a_{30}) +(1+\eta ^2)^2}{(1+2\eta ^2)+(1-\eta ^4)(a_{03}+a_{30})+(1-\eta ^2)^2} \\ \lambda '_{2,3}&=\frac{2\eta ^2(\rho _{23} \pm |\rho _{14}|)}{\rho _{11}+\eta ^2(\rho _{22}+\rho _{33})+\eta ^4 \rho _{44}} \\&=\frac{2\eta ^2 (a_{+-}+|a_{--}|)}{(1+2\eta ^2)+(1-\eta ^4)(a_{03}+a_{30})+(1-\eta ^2)^2} \end{aligned} \end{aligned}$$

(8.11)

Now, to implement the violation of the Bell-CHSH inequality in de Sitter we follow the following steps:

• 1. Step 1:

  After passing through the local filter in the new basis we have to satisfy the following necessary constraint condition:

  $$ c'(c')^{\dagger }>1 $$

• 2. Step 2:

  The above condition directly implies the following inequality:

  $$ 0>\eta ^4-\frac{(\rho _{23} + |\rho _{14}|)^2}{\rho _{44} (\rho _{22}+\rho _{33})}\eta ^2 + \frac{\rho _{11}}{\rho _{44}} $$

  where each of the quantities is appearing previously in the density matrix after passing through local filter.

• 3. Step 3:

  The above inequality can further be simplified to

  $$ (\rho _{23}+|\rho _{14}|)^4> 4\rho _{11}\rho _{44}(\rho _{22}+\rho _{33})^2 $$

• 4. Step 4:

  For real parameter $\eta$, if we substitute the entries of the local filter transformed density matrix in terms of the time dependent Bloch coefficients in the $(+,-,3)$ transformed basis, then the above inequality becomes:

  $$ (a_{+-}+|a_{--}|)^4 >(1-a_{33})^2[(1+a_{33})^2-(a_{03}+a_{30})^2] $$

• 5. Step 5:

  Now, our job is to explicitly verify this inequality for our open quantum system described by two entangled atoms. To serve this purpose we plot the following functions:

  $$\begin{aligned} {{\mathcal {J}}}_1(t)= & {} (a_{+-}(t)+|a_{--}(t)|)^4 \end{aligned}$$

  (8.12)
  $$\begin{aligned} {\mathcal {J}}_2(t)= & {} (1-a_{33}(t))^2 \nonumber \\&[(1+a_{33}(t))^2-(a_{03}(t)+a_{30}(t))^2] \end{aligned}$$

  (8.13)

  separately. In the plot we represent ${{\mathcal {J}}}_1(t)$ and ${\mathcal {J}}_2(t)$ with green and red color.

The plots shows the behaviour of the two functions appearing on either sides of the Bell’s inequality. The above plots clearly shows that the function represented by green colour is always greater than that of red colour i.e the condition ${{\mathcal {J}}}_1(t)>{{\mathcal {J}}}_2(t)$ is always satisfied irrespective of the parameter chosen to test the inequality. This establishes Bell-CHSH inequality violation and non locality in De-Sitter Space with the present two atomic OQS setup.

9 General discussion

As seen above, all the entanglement measures give almost similar behaviour with respect to different parameters. In this section we will analyze those plots together.

In Figs. 7a, 8a, 9a, 10a, 11, 12a, 13, 14, 15a, 16a we have plotted different entanglement measures of our two atomic OQS set up in the static patch of de-Sitter space with respect to Time (T) keeping all other parameters fixed. Here T is defined as the time difference, $T=\tau -\tau ^{'}$, which is very useful for further analysis. Here $\tau$ is the usual time scale and $\tau ^{'}$ is the time scale where the equilibrium boundary condition is imposed. These plots have been normalized with $T=10^4$. It is clearly observed that initially the value of normalized measure is almost zero implying that our two atomic OQS set up do not show any signature of quantum entanglement. This is consistent with our assumption that initially there is no correlation and the quantum states being represented by pure states only. As time goes on, the sub-system gets more and more entangled due to more correlation and at a late time scale, it saturates to unity, the maximum value of the measure. The late time scale behaviour is also consistent with the prediction from the present system i.e. as time goes on the system is represented by mixed quantum states due to getting more quantum correlation.

In Figs. 7b, 8b, 9b, 10b, 12b, 15b, 16b we have plotted different entanglement measures with respect to frequency($\omega$). These plots have been normalized with $\omega =5\times 10^{-2}$. It is can be seen that initially for a finite frequency scale, the value of the measures are almost constant showing entanglement due to the appearance of mixed quantum states. As frequency is further increased the value of the measures decrease implying reduction in entanglement and when it gets close to $\omega _0$ the value of entropy is almost zero indicating that the joint density matrix of our subsystem becomes separable (Figs. 17, 18).

In Figs. 7c, 8c, 9c, 10c, 12c, 15c, 16c we have plotted different entanglement measures with respect to Euclidean distance(L). These plots have been normalized with ${ L}=5000$. It is clearly observed that initially there are some fluctuations with increase in entropy for small length scale. The maximum value of the fluctuation implies the maximum value of the entangled entropy one can obtain once we want to analyse the behaviour with respect to the Euclidean distance L. This is because at small L value the quantum states are dominated by pure states and there is no corresponding quantum correlation. But as the distance between the two atoms increases, the subsystem gets entangled to a maximum value close to 1. With further increase in distance, there is no change in entropy implying that the correlation between the atoms in our OQS set up is maximum and the corresponding quantum state is dominated by mostly mixed states. It can be observed that the measures increase with increase in the distance between the two atoms, indicating rise in entanglement with rise in L. As the distance between the two atoms increases, the subsystem gets entangled to a maximum value.

In Figs. 7d, 8d, 9d, 10d, 12d, 15d, 16d we have plotted different entanglement measures with respect to the parameter k. These plots have been normalized with $k=5\times 10^4$. We know that the curvature of static patch of the de Sitter space can be expressed in terms of the parameter k as:

$$\begin{aligned} R=\frac{12}{\alpha ^2}=\frac{12}{k^2+r^2}\approx \frac{12}{k^2}~~\mathrm{for}~~ k=\sqrt{\alpha ^2-r^2}\approx \alpha>>r\end{aligned}$$

(9.1)

which implies we actually have considered both flat and static de Sitter space by varying the parameter k. It is clearly observed that initially there are fluctuations in entanglement measures with increase in the value of k. However, the maximum value the fluctuation implies the maximum value of the entangled entropy one can obtain once we want to analyse the behaviour with respect to the parameter k. This is because at small k or non-zero effect of curvature value the quantum states are dominated by mixed states and the corresponding quantum correlation is non zero. With further increase in the value of the parameter k, the entanglement measures remain constant and very very small implying almost getting no quantum correlation and described by the pure quantum states for the large values of k, which corresponds to the flat space time situation.

Table 1 A comparative study of various entanglement measures for our two atomic OQS setup

Table 2 Comparative study of one atomic and two atomic OQS setup

In Table 1 a comparative study of different entanglement measure is done for our two atomic OQS setup. Various entanglement measures calculated in this context are compared taking into account various parameters like rescaled time, Frequency, Euclidean distance, Inverse curvature. The best entanglement measure is found after studying the entanglement with respect to the respective parameter in the entire chosen range and then the conclusion is given. In Table 2 a comparative study of the one atomic and two atomic systems are done and the main highlighting differences have been noted. The reduced subsystem density matrix shows entanglement between the two atoms constituting the system whereas the one atomic subsystem in not entangled. It is also noted that for the case of one atomic case, the equilibrium temperature of the bath is exactly equal to the Gibbons Hawking temperature whereas for the two atomic case the equilibrium temperature is written in terms of Gibbons Hawking and Unruh temperature. Due to this difference in temperature, the curvature of the background spacetime appears to be different for the two atomic cases. The Wightman function has only one component ($G_{11}$) for the one atomic case as the correlation function is independent of the position of the atom, but for the two atomic case the Wightman function has four components given by $G_{11}$, $G_{12}$, $G_{21}$ and $G_{22}$ corresponding to the dependence of the bath correlation functions on the position of the atoms. Entanglement measure observed (if any) in case of one atomic system is mainly due to its interactions with the bath whereas for the two atomic case significant amount of entanglement measures are observed due to their mutual entanglement besides interaction with the bath.

10 Conclusion

To summarize, in this work, we have addressed the following issues to study the quantum entanglement phenomenon from two entangled atomic OQS set up:

• To begin with, we have started our discussion with two entangled atomic OQS set up. In this framework, the two entangled atoms mimic the role of Unruh-De-Witt detectors, conformally coupled to a thermal bath which is modelled by a massless scalar field in this specific problem. Apart from that, within this OQS set up, a non adiabatic Resonant Casimir Polder Interaction(RCPI) takes place between the Unruh-De-Witt detectors and the thermal bath. Most importantly this interaction is effected by the background De-Sitter space time in which the probe massless scalar field is fluctuating.

• Since we are only interested in the dynamics of the reduced two entangled atomic subsystem, we partially trace over the probe massless scalar field or thermal bath degrees of freedom. Consequently, without solving the total(system+bath+interaction) quantum liouville equation for the total density matrix, we actually solve the Gorini-Kossakowski-Sudarshan-Lindblad equation (Master equation) to explicitly know about the time evolution of the reduced subsystem density matrix. However, solving GSKL master equation with proper initial condition is an extremely complicated task, as it involves two non trivial components in the equation of motion. These are the effective Hamiltonian and Quantum dissipator or Lindbladian operator. Due to the complicated structure of both of them it is obvious that the analytical solution of the GSKL master equation is not possible for all type of OQS setup. For our two atomic entangled OQS setup, we represent the density matrix corresponding to each of the atom through bloch sphere representation. However, the reduced subsystem density matrix, which can be constructed by taking the tensor product of two atomic density matrices, cannot be parametrized by a bloch sphere. Instead of that we found the reduced subsystem density matrix is actually parametrized by three time dependent coefficients $a_{0i}(\tau ), a_{i0}(\tau )$ and $a_{ij}(\tau )$ $\forall i,j=1,2,3$. this implies that solving GSKL master equation for the present OQS set up is actually determining the time dependent behaviour of the above mentioned coefficients. We found that this leads to huge number of coupled differential equations of all these time dependent coefficients, which are not analytically solvable for given appropriate initial condition. To solve this problem next we transformed the basis from $\left\{ 1,2,3 \right\}$ to $\left\{ +,-,3\right\}$. In this new basis, we get simplified form of the linear differential equation which are less in number compared to the previous case. Also using the large time equilibrium behaviour of the reduced density matrix, which plays the role of initial condition in our problem, we have found the explicit analytical solution of $a_{ij} \ \forall \ \ i,j= +,-,3$.

• Using the analytical density matrix of the reduced subsystem, we further computed various measures of quantum entanglement viz. Von Neumann entropy, R${e'}$nyi entropy, Logarithmic negativity, Quantum discord, Entanglement of formation and Concurrence, which are commonly used in the context of quantum information theory these days. From the time dependent behaviour of all the measure of quantum entanglement we have found out almost the similar behaviour which states that for initial time $t=0$, the measure of quantum entanglement is 0 and as time goes on the subsystem gets more entangled and after a certain time it increases very slowly i.e. it almost saturates. Apart from these as we have obtained the similar feature both in the case of Von Neumann entropy and in the case of Quantum discord, this imply existence of long range quantum correlation at the late time scale, satisfying the necessary and sufficient condition for quantum entanglement. Similarly we have obtained the time dependent feature of logarithmic negativity, entanglement of formation and concurrence which strongly imply that at initial time $t=0$ our two atomic OQS setup do not show any signature of quantum entanglement as in all the cases the quantum measure is zero. As time goes on we have found out all of these measures significantly increase and at very late time scale it almost saturates. This is obviously a significant finding of our two atomic entangled OQS setup from which one can extract the existence of long range quantum correlation in late time scale, which is a very common topic of study in the context of quantum information theory.

• Last but not the least we have studied Bell-CHSH inequality⁸ violation from our present setup in de-Sitter space. Though these kind of violation of Bell-CHSH inequality in de-Sitter space is not very trivial. Most importantly, without introducing any axion in a more model independent way we have established the violation of Bell-CHSH inequality in de-sitter space, which is the necessary ingredient to study the non local effect in correlations in quantum mechanics.

The future prospects of this work is as follows.

• In this work, we have restricted our subsystem which is made up of two entangled atoms. One can further generalize the same problem with arbitrary number of atoms within the framework of OQS.

• In this work, we did our computation in the background static De-Sitter metric. However this framework can be implemented in any curved spacetime metric. One can even carry forward the calculations in other patches of De-Sitter space like the global and inflationary (planer) patch. It is expected to get significant modifications in the cosmological correlation functions. Within the framework of OQS and particularly for inflationary patch [54,55,56,57,58,59] of the De-Sitter space, one can exactly compare these results with the known cosmological correlation functions computed within the framework of closed quantum system. Such comparative analysis between the obtained cosmological correlation functions obtained from OQS and CQS will help us to know about the correct quantum mechanical picture of early universe cosmology. Additionally, by comparing this result with the data obtained from the observational probe for the early universe cosmology one may further rule out one of the possible quantum pictures mentioned here.

• For our work, we have computed the signature of quantum entanglement from various quantum information theoretic measures. However, we have not done the calculation for all possible measures. For completeness and also to conclude about the existence of long range quantum correlations at the late time scale one can further compute squashed entangled entropy [60], entanglement of distillation [61], relative entropy [62] etc. Additionally, one can also compute fisher information from the present OQS set up to know about the difference between the results obtained in the classical and quantum limiting situations.

• This kind of study can be also be used in parameter estimation theory and quantum metrology to estimate various parameters relevant to the theory but cannot be treated as any quantum mechanical observable. [63]

• Very recently the phenomena of quantum teleportation has been observed for the first time with qutrit states [64, 65], which is obviously an outstanding finding in the context of quantum information theory. The authors have succeeded in teleporting three-dimensional quantum states for the first time. High-dimensional teleportation could play an important role in future quantum computers. In this direction, our future plan is to study such possibilities from the present OQS set up in de Sitter space.

Appendices

Wightman function for probe massless scalar field in static de Sitter Space

In this section we compute the two atomic Wightman correlation function for a massless probe scalar field in static de Sitter spacer characterised by the following infinitesimal line element:

$$\begin{aligned}&ds^{2}=\left( 1-\frac{r^{2}}{\alpha ^{2}}\right) dt^{2} -\left( 1-\frac{r^{2}}{\alpha ^{2}}\right) ^{-1}dr^{2} \nonumber \\&\quad -r^{2}(d\theta ^{2}+\sin ^{2}\theta d\phi ^{2})~~~\mathrm{where}~~\alpha =\sqrt{\frac{3}{\Lambda }}>0. \end{aligned}$$

(A.1)

To compute the expression for the each of the entries of the two body Wightman function of the probe scalar field present in the external thermal bath we use the four dimensional static de Sitter geometry of our space-time. In this set of coordinate system in four dimension, the Klein-Gordon field equation for the massless conformally coupled external probe scalar field for the non-adiabatic environment can be expressed as:

$$\begin{aligned}&\left[ \frac{1}{\cosh ^3\left( \frac{t}{\alpha }\right) }\frac{\partial }{\partial t}\left( \cosh ^3\left( \frac{t}{\alpha }\right) \frac{\partial }{\partial t}\right) -\frac{1}{\alpha ^2\cosh ^2\left( \frac{t}{\alpha }\right) }\mathbf{L}^2\right] \nonumber \\&\quad \Phi (t,\chi ,\theta ,\phi )=0, \end{aligned}$$

(A.2)

where $\mathbf{L}^2$ is the Laplacian differential operator in the three dimensions characterised by the coordinate $(\chi ,\theta ,\phi )$ , which is explicitly defined as:

$$\begin{aligned} \mathbf{L }^2= & {} \frac{1}{\sin ^2\chi } \nonumber \\&\left[ \frac{\partial }{\partial \chi } \left( \sin ^2\chi \frac{\partial }{\partial \chi }\right) +\frac{1}{\sin \theta }\frac{\partial }{\partial \theta } \left( \sin \theta \frac{\partial }{\partial \theta }\right) \right. \nonumber \\&\quad \left. +\frac{1}{\sin ^2\theta }\frac{\partial ^2}{\partial \phi ^2}\right] , \end{aligned}$$

(A.3)

where we introduce a new coordinate $\chi$ which is related to the radial coordinate r as $r=\sin \chi$. The corresponding two body Wightman function between two space-time points for massless probe scalar field can be expressed as:

$$\begin{aligned} G(x,x^{'})= & {} \begin{pmatrix} ~ G^{11}(x,x')~~~ &{}~~~ G^{12}(x,x')~ \\ ~G^{21}(x,x')~~~ &{}~~~ G^{22}(x,x') ~ \end{pmatrix} \nonumber \\= & {} \begin{pmatrix} ~ \langle \Phi (\mathbf{x_{1}},\tau )\Phi (\mathbf{x_{1}},\tau ')\rangle ~~~ &{}~~~ \langle \Phi (\mathbf{x_{1}},\tau )\Phi (\mathbf{x_{2}},\tau ')\rangle ~ \\ ~\langle \Phi (\mathbf{x_{2}},\tau )\Phi (\mathbf{x_{1}},\tau ')\rangle ~~~ &{}~~~ \langle \Phi (\mathbf{x_{2}},\tau )\Phi (\mathbf{x_{2}},\tau ')\rangle ~ \end{pmatrix}~\nonumber \\= & {} \frac{1}{Z_\mathbf{Bath }}\begin{pmatrix} ~ \mathrm{Tr}\left[ \rho _\mathbf{Bath }(\tau -\tau ^{'}) \Phi (\mathbf{x_{1}},\tau )\Phi (\mathbf{x_{1}},\tau ')\right] ~ &{}~\mathrm{Tr}\left[ \rho _\mathbf{Bath }(\tau -\tau ^{'}) \Phi (\mathbf{x_{1}},\tau )\Phi (\mathbf{x_{2}},\tau ')\right] ~ \\ ~\mathrm{Tr}\left[ \rho _\mathbf{Bath }(\tau -\tau ^{'}) \Phi (\mathbf{x_{2}},\tau )\Phi (\mathbf{x_{1}},\tau ')\right] ~ &{}~ \mathrm{Tr}\left[ \rho _\mathbf{Bath }(\tau -\tau ^{'}) \Phi (\mathbf{x_{2}},\tau )\Phi (\mathbf{x_{2}},\tau ')\right] ~ \end{pmatrix}, \end{aligned}$$

(A.4)

where the Partition function for the bath is given by:

$$\begin{aligned} Z_\mathbf{Bath }=\mathrm{Tr}\left[ \rho _\mathbf{Bath }(\tau -\tau ^{'})\right] . \end{aligned}$$

(A.5)

Also, the components of the two atomic Wightman function can be expressed as⁹:

$$\begin{aligned} G^{11}(x,x')= & {} G^{22}(x,x') =\langle \Phi (\mathbf{x_{1}},\tau )\Phi (\mathbf{x_{1}},\tau ')\rangle \nonumber \\= & {} \langle \Phi (\mathbf{x_{2}},\tau )\Phi (\mathbf{x_{2}},\tau ')\rangle \nonumber \\= & {} - \frac{1}{16\pi ^{2}k^{2}}\frac{1}{\sinh ^2\left( \frac{\Delta \tau }{2k}-i\epsilon \right) }, \end{aligned}$$

(A.6)

$$\begin{aligned} G^{12}(x,x')= & {} G^{21}(x,x') = \langle \Phi (\mathbf{x_{1}},\tau )\Phi (\mathbf{x_{2}},\tau ')\rangle \nonumber \\= & {} \langle \Phi (\mathbf{x_{2}},\tau )\Phi (\mathbf{x_{1}},\tau ')\rangle \nonumber \\= & {} -\frac{1}{16\pi ^{2}k^{2}}\frac{1}{\left\{ \sinh [2](\frac{\Delta \tau }{2k}-i\epsilon )-\frac{r^{2}}{k^{2}}\sin [2](\frac{\Delta \theta }{2})\right\} }. \end{aligned}$$

(A.7)

Here we have introduced few parameters, which are defined as:

$$\begin{aligned} k= & {} \sqrt{g_{00}}\alpha = \sqrt{\alpha ^{2}-r^{2}}, \nonumber \\ \Delta \tau= & {} \tau -\tau '=\sqrt{g_{00}}(t-t')=k\left( \frac{t-t'}{\alpha }\right) . \end{aligned}$$

(A.8)

Gorini Kossakowski Sudarshan Lindblad (GSKL) ($C^{\alpha \beta }_{ij})$ matrix

In this appendix, we explicitly write down the entries of the Gorini Kossakowski Sudarshan Lindblad (GSKL) ($C^{\alpha \beta }_{ij})$ matrix which is appearing in the expression for the quantum dissipator or Lindbladian operator as given by the following expression:

$$\begin{aligned}&\mathcal {L}[\rho _\mathbf{System }(\tau )] \nonumber \\&\quad =\frac{1}{2}\sum _{i,j=\pm }^3 \sum _{\alpha \beta =1}^2 C_{ij}^{\alpha \beta } \nonumber \\&\qquad \left[ 2 \sigma _j^{\beta }\rho _\mathbf{System }(\tau ) \sigma _i^\alpha -\left\{ \sigma _i^\alpha \sigma _j^\beta ,\rho _\mathbf{System }(\tau )\right\} \right] , \end{aligned}$$

(B.1)

where we have written the expression in a transformed basis span by $(+,-,3)$ for two entangled OQS set up. The components of $C^{\alpha \beta }_{ij}$ matrix are very crucial to solve the time evolution equation of the reduced subsystem density matrix when we partially trace out the bath degrees of freedom from our two atomic entangled OQS set up. In general, Gorini Kossakowski Sudarshan Lindblad (GSKL) ($C^{\alpha \beta }_{ij})$ can be expressed as:

$$\begin{aligned}&C^{\alpha \beta }_{ij}={\tilde{A}}^{\alpha \beta }\delta _{ij}-i{\tilde{B}}^{\alpha \beta }\epsilon _{ijk}\delta _{3k}-{\tilde{A}}^{\alpha \beta }\delta _{3i}\delta _{3j}\nonumber \\&\quad \forall ~~i,j=+,-,3~~\mathrm{and}~~\forall \alpha ,\beta =1(\mathbf{Atom~1}),2(\mathbf{Atom~2}).\nonumber \\ \end{aligned}$$

(B.2)

In terms of explicit components the entries of the $C^{\alpha \beta }_{ij}$ matrix can be written as:

$$\begin{aligned}&C_{++}^{\alpha \beta }={\tilde{A}}^{\alpha \beta }, ~ C_{+-}^{\alpha \beta }=-i {\tilde{B}}^{\alpha \beta }, \nonumber \\&C_{+3}^{\alpha \beta }= C_{-3}^{\alpha \beta }=C_{3+}^{\alpha \beta }=C_{3-}^{\alpha \beta }=C_{33}^{\alpha \beta }=0, \nonumber \\&C_{-+}^{\alpha \beta }=i {\tilde{B}}^{\alpha \beta },~ C_{--}^{\alpha \beta }= {\tilde{A}}^{\alpha \beta } \end{aligned}$$

(B.3)

In the next two subsections, we will explicitly compute each of the components of ${\tilde{A}}^{\alpha \beta }$ and ${\tilde{B}}^{\alpha \beta }$ for the two atomic OQS in $(+,-,3)$ transformed basis.

Now we are interested in the following limit¹⁰:

$$\begin{aligned}&\mathrm{coth}(\pi k \omega _0)=0~\Longrightarrow ~k\omega _0=i\left( n+ \frac{1}{2}\right) \nonumber \\&\quad \mathrm{where}~n\in {\mathbb {Z}} \end{aligned}$$

(B.4)

Finally, substituting the explicit forms of ${\tilde{A}}^{\alpha \beta }$ and ${\tilde{B}}^{\alpha \beta }$ which we have derived in the previous sub section, we get the following expressions for the entries of the GSKL matrix:

$$\begin{aligned} C_{+3}^{\alpha \beta }= & {} C_{-3}^{\alpha \beta }=C_{3+}^{\alpha \beta }=C_{3-}^{\alpha \beta }=C_{33}^{\alpha \beta }=0 \end{aligned}$$

(B.5)

$$\begin{aligned} C_{-+}^{11}= & {} C_{-+}^{22}=-C_{+-}^{11}=-C_{+-}^{22}\nonumber \\= & {} i {\tilde{B}}^{11}=i{\tilde{B}}^{22}=\frac{\mu ^2\omega _0 i}{8\pi }, \end{aligned}$$

(B.6)

$$\begin{aligned} C_{-+}^{12}= & {} C_{-+}^{21}=-C_{+-}^{12}=-C_{+-}^{21}=i {\tilde{B}}^{12}\nonumber \\= & {} i {\tilde{B}}^{21}=\frac{\mu ^2\omega _0 i}{8\pi }~f\left( \omega _0,\frac{L}{2}\right) \end{aligned}$$

(B.7)

$$\begin{aligned} C_{++}^{11}= & {} C_{++}^{11}=C_{--}^{11}=C_{--}^{22}\nonumber \\= & {} {\tilde{A}}^{11}= {\tilde{A}}^{22}=\frac{\mu ^2\omega _0}{8\pi }\mathrm{coth}\left( \pi k\omega _0\right) , \end{aligned}$$

(B.8)

$$\begin{aligned} C_{++}^{12}= & {} C_{++}^{21}=C_{--}^{12}=C_{--}^{21}\nonumber \\= & {} {\tilde{A}}^{12}={\tilde{A}}^{21}=\frac{\mu ^2\omega _0}{8\pi }\mathrm{coth}\left( \pi k\omega _0\right) ~f\left( \omega _0,\frac{L}{2}\right) . \end{aligned}$$

(B.9)

Further, using the assumption, $k\omega _0=i\left( n+ \frac{1}{2}\right)$ where $n\in {\mathbb {Z}}$, we get the following simplified expressions for the entries of the GSKL matrix:

$$\begin{aligned} C_{+3}^{\alpha \beta }= & {} C_{-3}^{\alpha \beta }=C_{3+}^{\alpha \beta }=C_{3-}^{\alpha \beta }=C_{33}^{\alpha \beta }=0 \end{aligned}$$

(B.10)

$$\begin{aligned} C_{-+}^{11}= & {} C_{-+}^{22}=-C_{+-}^{11}=-C_{+-}^{22}=i {\tilde{B}}^{11}\nonumber \\= & {} i {\tilde{B}}^{22}=-\frac{\mu ^2}{8\pi k}\left( n+\frac{1}{2}\right) , \end{aligned}$$

(B.11)

$$\begin{aligned} C_{-+}^{12}= & {} C_{-+}^{21}=-C_{+-}^{12}=-C_{+-}^{21}=i {\tilde{B}}^{12} \nonumber \\= & {} i {\tilde{B}}^{21}=-\frac{\mu ^2}{8\pi k}\left( n+\frac{1}{2}\right) ~f\left( \omega _0,\frac{L}{2}\right) \end{aligned}$$

(B.12)

$$\begin{aligned} C_{++}^{11}= & {} C_{++}^{11}=C_{--}^{11}=C_{--}^{22}= {\tilde{A}}^{11}= {\tilde{A}}^{22}=0,\end{aligned}$$

(B.13)

$$\begin{aligned} C_{++}^{12}= & {} C_{++}^{21}=C_{--}^{12}=C_{--}^{21}= {\tilde{A}}^{12}={\tilde{A}}^{21}=0, \end{aligned}$$

(B.14)

where, the parameter

$$\begin{aligned} k=\sqrt{\frac{3}{\Lambda }-r^2}=\sqrt{\left( \frac{12}{R_\mathbf{DS }}\right) ^2-r^2}>0, \end{aligned}$$

(B.15)

for de Sitter space. Here we use the fact that, the curvature of static de Sitter space is given by, $R_\mathbf{DS }=\sqrt{48\Lambda }>0$ as $\Lambda >0$ for De Sitter.

Effective Hamiltonian ($H_{ij}^{\alpha \beta }$) matrix

In this appendix, we explicitly write down the entries of the Effective Hamiltonian ($H_{ij}^{\alpha \beta }$) matrix which is appearing in the expression for the Lamb Shift part of the Hamiltonian as given by the following expression:

$$\begin{aligned} H_\mathbf{Lamb Shift }=-\frac{i}{2}\sum ^{2}_{\alpha ,\beta =1}\sum ^{3}_{i,j=1}H_{ij}^{\alpha \beta } (n^{\alpha }_{i}.\sigma ^{\alpha }_{i})(n^{\beta }_{j}.\sigma ^{\beta }_{j}). \end{aligned}$$

(C.1)

where this Hamiltonian is constructed by partially tracing over bath degrees of freedom, which is in our problem a probe massless scalar field. Technically, the Effective Hamiltonian ($H_{ij}^{\alpha \beta }$) matrix represent here the strength of the spin chain interaction mentioned above. The general expression of $H_{ij}^{\alpha \beta }$ matrix is given by the following expression:

$$\begin{aligned} H^{ij}_{\alpha \beta }= & {} {\mathcal {A}}^{\alpha \beta } \delta _{ij}-i {\mathcal {B}}^{\alpha \beta } \epsilon _{ijk}\delta _{3k}-{\mathcal {A}}^{\alpha \beta } \delta _{3i} \delta _{3j}\nonumber \\ \forall ~~i,j= & {} +,-,3~~\mathrm{and}~~\forall \alpha ,\beta =1(\mathbf{Atom~1}),2(\mathbf{Atom~2}).\nonumber \\ \end{aligned}$$

(C.2)

In terms of explicit components the entries of the ${H}_{ij}^{\alpha \beta }$ matrix can be written as:

$$\begin{aligned} H_{++}^{\alpha \beta }= & {} H_{--}^{\alpha \beta }=A^{\alpha \beta }~~~\forall ~\alpha ,\beta =1,2, \end{aligned}$$

(C.3)

$$\begin{aligned} H_{+-}^{\alpha \beta }= & {} -i B^{\alpha \beta }~~~\forall ~\alpha ,\beta =1,2, \end{aligned}$$

(C.4)

$$\begin{aligned} H_{-+}^{\alpha \beta }= & {} i B^{\alpha \beta }~~~\forall ~\alpha ,\beta =1,2, \end{aligned}$$

(C.5)

$$\begin{aligned} H_{+3}^{\alpha \beta }= & {} H_{-3}^{\alpha \beta }=0~~~\forall ~\alpha ,\beta =1,2, \end{aligned}$$

(C.6)

$$\begin{aligned} H_{3j}^{\alpha \beta }= & {} 0~~~\forall ~j=+,-,3~~\mathrm{and}~~\forall ~\alpha ,\beta =1,2. \end{aligned}$$

(C.7)

Finally, substituting the explicit forms of ${A}^{\alpha \beta }$ and ${B}^{\alpha \beta }$ which we have derived in the previous sub section, we get the following expressions for the entries of the effective Hamiltonian matrix:

$$\begin{aligned} H_{ij}^{11}= & {} H_{ij}^{22}=\frac{\mu ^2P}{4\pi ^2 i}\nonumber \\&\left[ (\delta _{ij}-\delta _{3i}\delta _{3j})\Theta _1-i \epsilon _{ijk}\delta _{3k}\Theta _2 \right] ,\end{aligned}$$

(C.8)

$$\begin{aligned} H_{ij}^{12}= & {} H_{ij}^{21}=\frac{\mu ^2P}{4\pi ^2 i}\nonumber \\&\left[ (\delta _{ij}-\delta _{3i}\delta _{3j})\Theta _3-i \epsilon _{ijk}\delta _{3k}\Theta _4 \right] . \end{aligned}$$

(C.9)

Therefore, all the various entries of ${H}_{ij}^{\alpha \beta }$ can be written as:

$$\begin{aligned} H_{++}^{11}= & {} H_{++}^{22}=H_{--}^{11}=H_{--}^{22} \nonumber \\= & {} A^{11}=A^{22}=\frac{\mu ^2 P}{4\pi ^2 i}\Theta _1 \end{aligned}$$

(C.10)

$$\begin{aligned} H_{-+}^{11}= & {} H_{-+}^{22}=-H_{+-}^{11}=-H_{+-}^{22}\nonumber \\= & {} i B^{11}=i B^{22}=\frac{\mu ^2 P}{4\pi ^2 }\Theta _2 \end{aligned}$$

(C.11)

$$\begin{aligned} H_{++}^{12}= & {} H_{++}^{21}=H_{--}^{12}=H_{--}^{21}\nonumber \\= & {} A^{12}=A^{21}=\frac{\mu ^2 P}{4\pi ^2 i}\Theta _3 \end{aligned}$$

(C.12)

$$\begin{aligned} H_{-+}^{12}= & {} H_{-+}^{21}=-H_{+-}^{12}=-H_{+-}^{21}\nonumber \\= & {} i B^{12}=i B^{21}=\frac{\mu ^2 P}{4\pi ^2 }\Theta _4 \end{aligned}$$

(C.13)

The function $f\left( \omega ,\frac{L}{2}\right)$ has already been defined in the previous section. Here we introduce four integral functions $\Theta _i\forall i=1,2,3,4$, which is defined as:

$$\begin{aligned}&{{{\underline{\varvec{Integral I:}}}}}~~~ \Theta _{1}: \nonumber \\&\quad =\int ^{\infty }_{-\infty }d\omega ~\frac{\omega ^2}{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) } \end{aligned}$$

(C.14)

$$\begin{aligned}&{{{\underline{\varvec{Integral II:}}}}}~~~ \Theta _{2}:\nonumber \\&\quad =\int ^{\infty }_{-\infty }d\omega ~\frac{\omega \omega _0}{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) } \end{aligned}$$

(C.15)

$$\begin{aligned}&{{{\underline{\varvec{Integral III:}}}}}~~~ \Theta _{3}: \nonumber \\&\quad =\int ^{\infty }_{-\infty }d\omega ~\frac{\omega ^2 }{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) }f\left( \omega ,\frac{L}{2}\right) \nonumber \\ \end{aligned}$$

(C.16)

$$\begin{aligned}&{{{\underline{\varvec{Integral IV:}}}}}~~~ \Theta _{4}: \nonumber \\&\quad =\int ^{\infty }_{-\infty }d\omega ~\frac{\omega \omega _0}{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) }f\left( \omega ,\frac{L}{2}\right) .\nonumber \\ \end{aligned}$$

(C.17)

In the next section we explicitly compute the contributions from all of these integrals.

Calculation of useful integrals

In this section, we explicitly compute the analytical expression for the useful integrals $\Theta _i\forall i=1,2,3,4$, which are very useful to compute the expressions for the effective Hamiltonian matrix elements.

1.1 Integral I

In this subsection we explicitly compute the finite contribution from the following integral:

$$\begin{aligned} {{\Theta _{1}:=\int ^{\infty }_{-\infty }d\omega ~\frac{\omega ^2}{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) }}}~. \end{aligned}$$

(D.1)

In the limiting approximation $2\pi k\omega>>1$, one can further expand the integrand by taking large $k\omega$ approximation as:

$$\begin{aligned} {{\mathcal {V}}}(\omega _0,\omega ,k)&:=~\frac{\omega ^2}{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) } \nonumber \\&\xrightarrow [2\pi k\omega>>1]{}\frac{\omega ^2}{\left( \omega +\omega _0\right) \left( \omega -\omega _0\right) } :={\mathcal {V}}(\omega _0,\omega ). \end{aligned}$$

(D.2)

This implies that, after taking large $k\omega$ approximation the integrand of $\Theta _1$ becomes independent of the parameter k. Now, further using this approximation the integral $\Theta _1$ can be further simplified as:

$$\begin{aligned} \Theta _{2}\approx & {} \int ^{\infty }_{-\infty }d\omega ~{\mathcal {V}}(\omega _0,\omega ) \nonumber \\= & {} \underbrace{\int ^{0}_{-\infty }d\omega ~{\mathcal {V}}(\omega _0,\omega )}_{{\equiv ~ {\mathcal {D}}_1(\omega _0)}} ~+~\underbrace{\int ^{\infty }_{0}d\omega ~{\mathcal {V}}(\omega _0,\omega )}_{{\equiv ~ {{\mathcal {D}}}_2(\omega _0)}}, \end{aligned}$$

(D.3)

where we have written the integrals into two parts,represented by ${{{\mathcal {D}}}_1(\omega _0)}$ and ${{{\mathcal {D}}}_2(\omega _0)}$. Now, here we see that in the $2\pi k>>1$ limit we get:

$$\begin{aligned} {{{\mathcal {D}}}_1(\omega _0)}= & {} \int ^{0}_{-\infty }d\omega ~{\mathcal {V}}(\omega _0,\omega ) \nonumber \\= & {} -\int ^{\infty }_{0}d\omega ~{\mathcal {V}}(\omega _0,\omega )=-{{{\mathcal {D}}}_2(\omega _0)}. \end{aligned}$$

(D.4)

Now, here ${{{\mathcal {D}}}_1(\omega _0)}$ and ${{{\mathcal {D}}}_2(\omega _0)}$ gives divergent contributions in the frequency range, $-\infty<\omega <0$ and $0<\omega <\infty$. To get the finite regularised contributions from these integrals we introduce a cut-off regulator $\omega _c$, by following Bethe regularisation procedure. After introducing this cut-off we get the following result:

$$\begin{aligned} {{{\mathcal {D}}}_1(\omega _0,\omega _c)}= & {} \int ^{0}_{-\omega _c}d \omega ~{\mathcal {V}}(\omega _0,\omega )=\int ^{\omega _c}_{0}d \omega ~{\mathcal {V}}(\omega _0,\omega ) \nonumber \\= & {} {{\mathcal {D}}_2(\omega _0,\omega _c)}=\frac{1}{2}\left[ \omega _c- \omega _0\tanh ^{-1}\left( \frac{\omega _c}{\omega _0}\right) \right] .\nonumber \\ \end{aligned}$$

(D.5)

Consequently, we get the following regularised expression for the integral $\Theta _1$, as given by:

$$\begin{aligned} \Theta _{1}= & {} {{{\mathcal {D}}}_1(\omega _0,\omega _c)}+ {{\mathcal {D}}_2(\omega _0,\omega _c)} \nonumber \\= & {} \left[ \omega _c-\omega _0\tanh ^{-1}\left( \frac{\omega _c}{\omega _0}\right) \right] . \end{aligned}$$

(D.6)

Now, if we further use the approximation that the cut-off is small compared to $\omega _0$ i.e. $\omega _c<<\omega _0$, then we get¹¹:

$$\begin{aligned} {{\boxed {\boxed {\mathbf{Integral ~I:}~~~\Theta _{2}={{\mathcal {D}}_1(\omega _0,\omega _c)}+{{\mathcal {D}}_2(\omega _0,\omega _c)}=\left[ \omega _c-\omega _0\left( \frac{\omega _c}{\omega _0}\right) \right] \approx 0}}}}~.\nonumber \\ \end{aligned}$$

(D.9)

1.2 Integral II

In this subsection we explicitly compute the finite contribution from the following integral:

$$\begin{aligned} {\Theta _{2}:=\int ^{\infty }_{-\infty }d\omega ~\frac{\omega \omega _0}{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) }}~. \end{aligned}$$

(D.10)

It is important to note that in the limit, $2\pi k\omega>>1$, one can further expand the integrand given by the following approximation as:

$$\begin{aligned} \begin{aligned} {\mathcal {M}}(\omega _0,\omega ,k)&:=~\frac{\omega _0~\omega }{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) } \\&\xrightarrow [2\pi k\omega>>1]{}\frac{\omega _0~\omega }{\left( \omega +\omega _0\right) \left( \omega -\omega _0\right) }:={\mathcal {M}}(\omega _0,\omega ). \end{aligned} \end{aligned}$$

(D.11)

This implies that, after taking limit $2\pi k\omega>>1$ the integrand of $\Theta _2$ becomes independent of the parameter k. Now, further using this approximation the integral $\Theta _1$ can be expressed as:

$$\begin{aligned} \Theta _{1}\approx & {} \int ^{\infty }_{-\infty }d\omega ~{\mathcal {M}}(\omega _0,\omega ) \nonumber \\= & {} \underbrace{\int ^{0}_{-\infty }d\omega ~{\mathcal {M}}(\omega _0,\omega )}_{{\equiv ~ {\mathcal {N}}_1(\omega _0)}} \nonumber \\&+ \underbrace{\int ^{\infty }_{0}d\omega ~{\mathcal {M}}(\omega _0,\omega )}_{{\equiv ~ {{\mathcal {N}}}_2(\omega _0)}}, \end{aligned}$$

(D.7)

where we have written the integrals into two parts, indicated by ${{{\mathcal {N}}}_1(\omega _0)}$ and ${{{\mathcal {N}}}_2(\omega _0)}$. Now, in the limit $2\pi k\omega>>1$, we get:

$$\begin{aligned} {{{\mathcal {N}}}_1(\omega _0)}= & {} \int ^{0}_{-\infty }d\omega ~{\mathcal {M}}(\omega _0,\omega ) \nonumber \\= & {} -\int ^{\infty }_{0}d\omega ~{\mathcal {M}}(\omega _0,\omega )=-{{{\mathcal {N}}}_2(\omega _0)}~. \end{aligned}$$

(D.12)

Now, here ${{{\mathcal {N}}}_1(\omega _0)}$ and ${{{\mathcal {N}}}_2(\omega _0)}$ gives divergent contributions in the frequency range, $-\infty<\omega <0$ and $0<\omega <\infty$. To get the finite contributions from these integrals we introduce a cut-off regulator $\omega _c$, by following Bethe regularisation procedure. After introducing this cut-off we get the following finite contribution:

$$\begin{aligned} {{{\mathcal {N}}}_1(\omega _0,\omega _c)}= & {} \int ^{0}_{-\omega _c}d\omega ~{\mathcal {M}}(\omega _0,\omega ) \nonumber \\= & {} -\int ^{\omega _c}_{0}d\omega ~{\mathcal {M}}(\omega _0,\omega ) \nonumber \\= & {} -{{\mathcal {N}}_2(\omega _0,\omega _c)} \nonumber \\= & {} -\frac{\omega _0}{2}\ln \left[ 1-\left( \frac{\omega _c}{\omega _0}\right) ^2\right] . \end{aligned}$$

(D.13)

Consequently, we get the following expression for the integral $\Theta _2$, as given by:

$$\begin{aligned}&{{\boxed {\boxed {\mathbf{Integral ~II:}~\Theta _{1}={{\mathcal {U}}_1(\omega _0,\omega _c)}+{{\mathcal {U}}_2(\omega _0,\omega _c)}=\frac{\omega _0}{2}\ln \left[ 1-\left( \frac{\omega _c}{\omega _0}\right) ^2\right] -\frac{\omega _0}{2}\ln \left[ 1-\left( \frac{\omega _c}{\omega _0}\right) ^2\right] =0}}}}. \end{aligned}$$

(D.14)

1.3 Integral III

In this subsection we explicitly compute the finite contribution from the following integral:

$$\begin{aligned} \begin{aligned} \Theta _{3}&:=\int ^{\infty }_{-\infty }d\omega ~\frac{\omega ^2}{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) }\\&\quad f\left( \omega ,\frac{L}{2}\right) ~, \end{aligned} \end{aligned}$$

(D.15)

where, we define the function $f\left( \omega ,\frac{L}{2}\right)$ given by the following expression:

$$\begin{aligned} {{f\left( \omega ,\frac{L}{2}\right) =\frac{1}{L\omega \sqrt{1+\left( \frac{L}{2k}\right) ^2}} \sin \left( 2k\omega \sinh ^{-1}\left( \frac{L}{2k}\right) \right) }}~.\nonumber \\ \end{aligned}$$

(D.16)

In the limit $2\pi k>>1$, one can further expand the integrand as:

$$\begin{aligned} \begin{aligned} {\mathcal {Z}}(\omega _0,\omega ,k)&:=~\frac{\omega ^2~f\left( \omega ,\frac{L}{2}\right) }{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) } \\&\xrightarrow [2\pi k\omega \>>1]{}\frac{\omega ^2~f\left( \omega ,\frac{L}{2}\right) }{\left( \omega +\omega _0\right) \left( \omega -\omega _0\right) }:={{\mathcal {Z}}(\omega _0,\omega ,k)}. \end{aligned} \end{aligned}$$

(D.17)

This implies that, in the limit $2\pi k>>1$ the integrand of $\Theta _3$ is not independent of the parameter k. Now, further using this approximation the integral $\Theta _3$ can be further simplified as:

$$\begin{aligned} \Theta _{3}\approx & {} \int ^{\infty }_{-\infty }d\omega ~{{\mathcal {Z}}(\omega _0,\omega ,k)} \nonumber \\= & {} \underbrace{\int ^{0}_{-\infty }d\omega ~\widetilde{{\mathcal {Z}}(\omega _0,\omega ,k)}}_{{\equiv ~ {\mathcal {R}}^{l}_1(\omega _0,k)}} \nonumber \\&+ \underbrace{\int ^{\infty }_{0}d\omega ~\widetilde{{\mathcal {Z}}(\omega _0,\omega ,k)}}_{{\equiv ~ {{\mathcal {R}}}^{l}_2(\omega _0,k)}}, \end{aligned}$$

(D.18)

where we have written the integrals into two parts, indicated by ${{{\mathcal {R}}}^{l}_1(\omega _0,k)}$ and ${{{\mathcal {R}}}^{l}_2(\omega _0,k)}$. Now, here in the limit $2\pi k>>1$ we get:

$$\begin{aligned} {{{\mathcal {R}}}^{l}_1(\omega _0,k)}= & {} \int ^{0}_{-\infty }d\omega ~{{{\mathcal {Z}}}(\omega _0,\omega ,k)} \nonumber \\= & {} \int ^{\infty }_{0}d\omega ~{{{\mathcal {Z}}}(\omega _0,\omega ,k)}={{{\mathcal {R}}}^{l}_2(\omega _0,k)}\nonumber \\= & {} \frac{\pi }{2L\sqrt{1+\left( \frac{L}{2k}\right) ^2}}\cos (2k\omega _0 \sinh ^{-1}\left( \frac{L}{2k}\right) )~.\nonumber \\ \end{aligned}$$

(D.19)

Consequently, we get the following expression for the integral $\Theta _3$, given by:

$$\begin{aligned} {{\boxed {\boxed {\mathbf{Integral ~III:}~~~\Theta _{3}={{\mathcal {R}}^{l}_1(\omega _0,k)}+{{\mathcal {R}}^{l}_2(\omega _0,k)}=\frac{\pi }{L\sqrt{1+\left( \frac{L}{2k}\right) ^2}}\cos (2k\omega _0 \sinh ^{-1}\left( \frac{L}{2k}\right) )}}}}~. \end{aligned}$$

(D.20)

Further, substituting $k\omega _0=i\left( n+\frac{1}{2}\right) ~\forall ~n\in {\mathbb {Z}}$, we get the following simplified expression for the integral $\Theta _3$, given by:

$$\begin{aligned} {{\boxed {\boxed {\mathbf{Integral ~III:}~~~\Theta _{3}={{\mathcal {R}}_1(\omega _0,k)}+{{\mathcal {R}}_2(\omega _0,k)}=\frac{\pi }{L\sqrt{1+\left( \frac{L}{2k}\right) ^2}}\cosh (\left( 2n+1\right) \sinh ^{-1}\left( \frac{L}{2k}\right) )}}}}~. \end{aligned}$$

(D.21)

1.4 Integral IV

In this subsection we explicitly compute the finite contribution from the following integral:

$$\begin{aligned} \begin{aligned} \Theta _{4}&:=\int ^{\infty }_{-\infty }d\omega ~\frac{\omega \omega _0}{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) }\\&\quad f\left( \omega ,\frac{L}{2}\right) ~, \end{aligned} \end{aligned}$$

(D.22)

where, we define the function $f\left( \omega ,\frac{L}{2}\right)$ given by the following expression:

$$\begin{aligned} \begin{aligned} f\left( \omega ,\frac{L}{2}\right)&=\frac{1}{L\omega \sqrt{1+\left( \frac{L}{2k}\right) ^2}} \\&\quad \sin \left( 2k\omega \sinh ^{-1} \left( \frac{L}{2k}\right) \right) ~. \end{aligned} \end{aligned}$$

(D.23)

In the limit $2\pi k>>1$, one can further expand the integrand as:

$$\begin{aligned} \begin{aligned} {\mathcal {S}}(\omega _0,\omega ,k)&:=~\frac{\omega \omega _0~f\left( \omega ,\frac{L}{2}\right) }{\left( 1-e^{-2\pi k\omega }\right) \left( \omega +\omega _0\right) \left( \omega -\omega _0\right) } \\&\xrightarrow [2\pi k\omega \>>1]{}\frac{\omega \omega _0~f\left( \omega ,\frac{L}{2}\right) }{\left( \omega +\omega _0\right) \left( \omega -\omega _0\right) }:={{\mathcal {S}}(\omega _0,\omega ,k)}. \end{aligned} \end{aligned}$$

(D.24)

This implies that, in the limit $2\pi k>>1$ the integrand of $\Theta _3$ is not independent of the parameter k. Now, further using this approximation the integral $\Theta _4$ can be further simplified as:

$$\begin{aligned} \Theta _{4}\approx & {} \int ^{\infty }_{-\infty }d\omega ~{{\mathcal {S}}(\omega _0,\omega ,k)} \nonumber \\= & {} \underbrace{\int ^{0}_{-\infty }d\omega ~{{\mathcal {S}}(\omega _0,\omega ,k)}}_{{\equiv ~ {\mathcal {L}}_1(\omega _0,k)}} \nonumber \\&+ \underbrace{\int ^{\infty }_{0}d\omega ~{{\mathcal {S}}(\omega _0,\omega ,k)}}_{{\equiv ~ {{\mathcal {L}}}_2(\omega _0,k)}}, \end{aligned}$$

(D.25)

where we have written the integrals into two parts, indicated by ${{{\mathcal {L}}}_1(\omega _0,k)}$ and ${{{\mathcal {L}}}_2(\omega _0,k)}$. Now, here in the limit $2\pi k>>1$ we get:

$$\begin{aligned} {{{\mathcal {L}}}_1(\omega _0,k)}= & {} \int ^{0}_{-\infty }d\omega ~{{{\mathcal {S}}}(\omega _0,\omega ,k)} \nonumber \\= & {} -\int ^{\infty }_{0}d\omega ~{{{\mathcal {S}}}(\omega _0,\omega ,k)}=-{{{\mathcal {L}}}_2(\omega _0,k)}\nonumber \\= & {} -\frac{i\sqrt{\pi }}{2L\sqrt{1+\left( \frac{L}{2k}\right) ^2}}\underbrace{G_{1,3}^{2,1}\left( -\omega ^2_0 k^2 \sinh ^{-1}\left( \frac{L}{2 k}\right) ^2| \begin{array}{c} \frac{1}{2} \\ \frac{1}{2},\frac{1}{2},0 \\ \end{array} \right) }_{\mathbf{Meijer~G~function}}~.\nonumber \\ \end{aligned}$$

(D.26)

Consequently, we get the following expression for the integral $\Theta _4$, as given by:

$$\begin{aligned} {{\boxed {\boxed {\mathbf{Integral ~IV:}~~~\Theta _{4}={{\mathcal {L}}_1(\omega _0,k)}+{{{\mathcal {L}}}_2(\omega _0,k)}=0}}}}~. \end{aligned}$$

(D.27)